Cost benefit analysis of the federal interventions that aim to reduce the harms of opioid crisis

Final Project | EPIB676

Author

Mariam El Sheikh

Published

21/04/2023 T11:39

Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/03b_fun_calib_ntd.R"))
source(here("02_scripts/06_interventions_models.R"))
theme_set(theme_minimal())

1 Background

—What is the decision problem your model aims to inform? What is already known, and what are the key questions you hope to address? Why is this important?—–

Opioid overdose deaths and related harms epidemic began in Canada in 2016 and is still ongoing. As a result, through the Substance Use and Addictions Program (SUAP) and other related programs, the government is funding projects that support people who use drugs and aim to mitigate the harms and deaths resulting from opioid overdose epidemic. For this project, I aim to evaluate the future economic impacts of two of these projects and their effect on reducing opioid-related deaths and harms.

This project aims to evaluate the future economic impact of these interventions by conducting a cost-benefit analysis to examine the difference in healthcare costs between the interventions and their projected effect on reducing opioid-related deaths and other important health outcomes

2 Methods

This project is a team project, I worked on the modeling part but many of the decisions and the work upon which the modeling was built was done by other teams members (DP, MS, IW). The model was designed and conceptualized by DP. A systematic review was conducted by DP, MS and IW to extract the effects of the interventions from the literature. Transition probabilities, and healthcare costs were extracted by DP from either the literature or were estimated based on scientific assumptions in consultation with experts working in the field and/or people with lived experiences. Many of the decisions made for the modelling part were made in consultation with DP.

Fig.1: A simple chematic of the Markov model representing pathways of opioid use in Canada

An open cohort Markov model was used (fig. 1) that consists of 31 states representing the pathways of opioid use in Canada over 15 years (2015 to 2029) with monthly cycles for 15+ Canadian population (designed and conceptualized by DP). It is an open cohort model that takes into account the changes in population over time by adding people to the first state (Pain free, no use) every cycle. The yearly increase was assumed to be constant throughout the year (meaning that each increase was divided by 12, as for years in the future, projections of population size in Canada were used and proportion of those 15+ was estimated based on previous data, then a similar approach was used to estimate the monthly increase in the cohort).

Table 1: List of variables corresponding to the states in the model

Code
flextable(ini_states_tbl_org %>% 
            select(`Description `, `Variable `)) |>
  autofit()

Description 

Variable 

Pain free, no use 

BN_PN

Acute pain, no use 

BN_ACUTE

Chronic (non-cancer) pain, no use 

BN_CHRONIC

Cancer, no use 

BN_CANCER

Other (i.e., cough, diarrhea, sickle cell disease), no use 

BN_OTHER

Palliative, no use 

BN_PALLIATIVE

Acute pain, Rx use 

BPO_ACUTE

Chronic (non-cancer) pain, Rx use 

BPO_CHRONIC

Cancer, Rx use 

BPO_CANCER

Other (i.e., cough, diarrhea, sickle cell pain), Rx use 

BPO_OTHER

Palliative, Rx use 

BPO_PALLIATIVE

Rx opioid misuse 

BPO_MISUSE

Illicit opioid use 

BI_ILLICIT

Detox / Withdrawal Management 

BS_DETOX

OAT initiation 

BS_OAT_INI

OAT maintenance 

BS_OAT_MAINT

OAT / Safe Supply 

BS_OAT_SS

Safe Supply 

BS_SS

Overdoes - illicit

BO_OD_ILLICIT

Overdose - misuse

BO_OD_RX

Moderate brain injury 

BO_MOD_BI

Severe brain injury 

BO_SEVERE_BI

Severe brain injury Out

BO_SEVERE_BI_OUT

Opioid-related death

BO_OD_DEATH

Death 

BO_DEATH

R: Illicit opioid use 

BR_ILLICIT

R: OAT initiation 

BR_OAT_INI

R: OAT maintenance 

BR_OAT_MAINT

R: OAT / Safe Supply 

BR_OAT_SS

R: Safe Supply 

BR_SS

R: Illicit opioid overdose

BR_OD_ILLICIT

Characteristics of this model:

  • One of the unique aspects of this model is that it distinguishes between chronic noncancer pain, and cancer pain and it includes palliative care, which are important pathways that we believe would be part of the spillover effects of at least one of the interventions we are examining (prescription guidelines).
  • We included brain injuries as part of the OD sequelea. We didn’t distinguish between mild brain injury and no brain injury, both transition from OD back into the general population, however we distinguished between moderate and severe brain injury. For severe brain injury as a result of an opioid-related OD we assumed that these individuals won’t be part of the model anymore, they can either transition to death or transition to an “out” state where they can stay there until they die (we created the “out” state so that the high brain injury healthcare costs wouldn’t be incurred every cycle for those who didn’t die). For the moderate brain injury state, we created another set of repeat states for opioid illicit use, substance use treatment that they are different from the general population states since we believed that the transitions would be different for those who experienced moderate brain injury.
  • It distinguishes between opioid-related OD deaths and deaths from any other causes.
  • It also distinguishes between Rx opioid misuse and opioid illicit use, and overdose as a result of an Rx opioid and opioids from illicit drug market.
  • It accounts for fentanyl contamination in the illicit drug market by including a multiplier of (1.005)^(t-36) (t-36 is number of month) to the probabilities of transitioning from opioid illicit use to illicit use OD and transitioning from illicit use OD to OD death, from 2018 to 2020. Then beyond 2020 it remained constant.

2.1 Alternatives

This model was run for 5 scenarios in regards to the previously mentioned interventions: 1) No interventions, 2) Naloxone only, 3) Prescription guidelines only, 4) Safer Supply only, and 5) all 3 interventions. The effects of each of the interventions were estimated by other team members who conducted a systematic review of the literature (DP, MS, IW). Transition probabilities were changed for the naloxone and prescription guidelines interventions (as a result the 1-rest transition probability will be changed).

For the naloxone intervention the following transition probabilities were changed starting from January 2017:

  • from illicit opioid OD to OD death was reduced by 5% (1-rest: illicit opioid OD to illicit opioid use), and
  • from rx opioid overdose to OD deaths was reduced by 3% (1-rest: Rx opioid OD to Rx opioid misuse).

For the prescription guidelines the following transition probabilities were changed starting from May 2017:

  • pain free to acute pain, rx use was reduced by 43% (1-rest: pain free stay in pain free),
  • cancer no use to cancer Rx use was reduced by 55% (1-rest: cancer no use stay in cancer no use),
  • acute pain Rx use to chronic pain, Rx use was reduced by 22% (1-rest: acute pain Rx use to pain free), and
  • chronic pain Rx use to chronic pain no use was increased by 4% (1-rest: chronic pain Rx use stay in chronic pain Rx use).

For safer supply to account for the effect of the pilot programs introduced in 2016, 1000 people were moved from illicit drug use state to OAT/Safer supply state in January 2016. Then to account for compassionate prescribing due to the COVID-19 panademic 10,000 people were moved from illicit drug use state to OAT/Safer supply state in April 2020. Before January 2016 all transitioning into and out of safer supply states (Safer Supply, OAT/Safer Supply, Repeat OAT/Safer Supply and Repeat Safer Supply). Starting from January 2016, all transitioning out of OAT/Safer Supply (and staying in) were changed to their corresponding non-zero values. Until in 2023, when the intervention is supposed to be implemented, all transitioning probabilities into and out of safer supply states were changed to their non-zero values.
For the cost-benefit analysis the effects of these interventions individually and all together was on interest, therefore 4 different comparisons were made: 1) Naloxone vs No Interventions, 2) Safer Supply vs No Interventions, 3) Prescription Guidelines vs No Interventions and 4) All Interventions vs No Interventions.

2.2 Outcomes

Primary Outcomes:

  • CBA: Opioid-related costs from the healthcare perspectives for all states were used in this analysis. These include costs for hospitalization, ED visits, physician billing, paramedic services, drugs, and other healthcare system related costs. These were identified from publicly available data, as well as published systematic reviews, studies, and reports. All costs were discounted using 0.03 discount factor annually.
  • Epidemiological measures: cumulative overall deaths and OD deaths.

Secondary Outcomes - more epidemiological measures:

  • Weighted yearly prevalence of 1) severe brain injuries due to opioid overdose, 2) moderate brain injuries due to opioid overdose, 3) substance use treatment and 4) prescription opioid use. In addition to trends in monthly prevalence counts for 4 outcomes.

For the primary outcomes changes and percent changes were calculated for the 4 different comparison groups while the secondary outcomes were calculated for each of th 5 scenarios.

2.3 Calibration and Uncertainty analysis

After designing the model structure, initial populations, and transitions, we calibrated the model using data from 2015 to 2021. Our calibration targets were total deaths (excluding opioid-related overdose deaths), the total number of people in Opioid Agonist Treatment (OAT), and proportions of prescription opioid use. To obtain this data, we used various sources, including Statistics Canada, the Canadian Institute for Health Information’s report on Opioid prescribing in Canada, and scientific literature (DP).
Even though we had opioid-related OD deaths numbers, we didn’t use them as targets because they are directly affected by different interventions that has been implemented so we couldn’t calibrate the no intervention parameters to them. However, we used them to assess if the estimates of opioid-related OD deaths from the model of the no intervention scenario are sound (they should be more than the targets).

We selected parameters to calibrate based on a one-way sensitivity analysis, where we changed one transition probability at a time, using lower and upper ranges, and observed how this impacted predicted outcomes (overall deaths and OD deaths). The one-way sensitivity analysis was conducted twice, 1) accounting for the fentanyl contamination (Scenario 1) and 2) without accounting for fentanyl contamination (Scenario 2). We included transition probabilities that resulted in a notable change in predicted outcomes (defined as 2.5%) in the set of parameters to be calibrated, in addition to all transition probabilities to death or opioid-related overdose deaths (a total of 38 parameters for scenario 1 and 39 paramters for scenario 2). We then conducted a Maximum A Posteriori (MAP) estimation using Nelder Mead, a sampling method that aims to find the global optimum twice once with the model for both scenarios. This allowed us to select the parameter set that best approximated the targets and provided us with information about changes in the prior information for some of the parameters. Overall, we updated the following baseline transition probabilities:

  • Pain free to Rx opioids misuse changed from 0.00032511 to 0.00220011
  • Chronic pain, rx use to rx OD from 0.00097898 to 0.00160398
  • Detox to illicit opioid use from 0.1818112 to 0.22181118
  • Detox to OAT initiation from 0.80000000 to 0.82000000
  • Palliative care, rx use to death from 0.37500000 to 0.38250000
  • OAT maintenance stay in OAT maintenance from 0.70000000 to 0.71000000

Then we used the Sample Importance Resampling (SIR) Bayesian calibration method after updating the prior information based on the MAP results, and we ran the model using the sample sets (that were randomly sampled from the priors defined for the parameters), estimated likelihood for each set, and weighted the sets by likelihood to approximate posteriors.

Probabilistic sensitivity analysis (PSA) was conducted on primary outcomes to account for the uncertainty in our result. The SIR generated posteriors for the parameters that were calibrated were used as priors and for other transition probabilities, a uniform distribution was used with ranges either from the literature and if not available, +/-5% were used as ranges. From these priors probabilistic samples were randomly samples and the analysis was simulated 10 000 times.

3 Results

What are your main findings? Should estimate some outcome(s) for >1 alternative, and potentially conduct incremental analysis (compute ICERS) depending on your framework. Report uncertainty analysis findings to give some idea of the robustness of your findings

3.1 Comparison between SIR sample, uncalibrated sample, uncalibrated point estimates and calibration targets

Fig.2a: Prevalence of prescription opioid use (an interactive version of this plot can be found in Appendix 4)

Fig.2b: Total number of overall deaths (excluding opioid-related OD deaths) (an interactive version of this plot can be found in Appendix 4)

Fig.2c: Yearly counts of people in OAT (an interactive version of this plot can be found in Appendix 4)

Fig.2d: Total number of opioid-related OD deaths (an interactive version of this plot can be found in Appendix 4)

The model approximated the calibration targets well: even though the distribution of prevalence of opioid prescription did not include the target values, the difference between the values was negligible. The distribution for overall deaths included the target values and the distribution for opioid-related deaths was larger than the target values (what we wanted the model to result in). However the distribution of OAT counts doesn’t include the target values and that might be either because these values aren’t accurate enough to represent the values estimated by the model or that the calibration process didn’t work as well. Overall, the estimates from the model seem to be approximating values found in the literature.

3.2 Results for the primary outcomes (PSA)

Code
m_outcomes_psa <- readRDS(here("01_data/m_outcomes_psa.RDS"))
point_est <- readRDS(here("01_data/point_estiamte.RDS"))
total_death_oddeaths <- readRDS(here("01_data/total_deaths_oddeaths.RDS"))


mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)

mean_cost_diff_nalox <- mean(m_outcomes_psa[, "cost_diff_nalox"])
ci_cost_diff_nalox <- quantile(m_outcomes_psa[, "cost_diff_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_ss <- mean(m_outcomes_psa[, "cost_diff_ss"])
ci_cost_diff_ss <- quantile(m_outcomes_psa[, "cost_diff_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_pg <- mean(m_outcomes_psa[, "cost_diff_pg"])
ci_cost_diff_pg <- quantile(m_outcomes_psa[, "cost_diff_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_all <- mean(m_outcomes_psa[, "cost_diff_all"])
ci_cost_diff_all <- quantile(m_outcomes_psa[, "cost_diff_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_cost_diff <- c(paste0(round(ci_cost_diff_nalox[3]/1000000, 0), " [",
                           round(ci_cost_diff_nalox[1]/1000000, 0),", ",
                           round(ci_cost_diff_nalox[2]/1000000, 0), "]"),
                    paste0(round(ci_cost_diff_ss[3]/1000000, 0), " [",
                           round(ci_cost_diff_ss[1]/1000000, 0),", ",
                           round(ci_cost_diff_ss[2]/1000000, 0), "]"),
                    paste0(round(ci_cost_diff_pg[3]/1000000, 0), " [",
                           round(ci_cost_diff_pg[1]/1000000, 0),", ",
                           round(ci_cost_diff_pg[2]/1000000, 0), "]"),
                    paste0(round(ci_cost_diff_all[3]/1000000, 0), " [",
                           round(ci_cost_diff_all[1]/1000000, 0),", ",
                           round(ci_cost_diff_all[2]/1000000, 0), "]"))


mean_death_diff_nalox <- mean(m_outcomes_psa[, "death_diff_nalox"])
ci_death_diff_nalox <- quantile(m_outcomes_psa[, "death_diff_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_ss <- mean(m_outcomes_psa[, "death_diff_ss"])
ci_death_diff_ss <- quantile(m_outcomes_psa[, "death_diff_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_pg <- mean(m_outcomes_psa[, "death_diff_pg"])
ci_death_diff_pg <- quantile(m_outcomes_psa[, "death_diff_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_all <- mean(m_outcomes_psa[, "death_diff_all"])
ci_death_diff_all <- quantile(m_outcomes_psa[, "death_diff_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_death_diff <- c(paste0(round(ci_death_diff_nalox[3], 0), " [",
                           round(ci_death_diff_nalox[1], 0),", ",
                           round(ci_death_diff_nalox[2], 0), "]"),
                    paste0(round(ci_death_diff_ss[3], 0), " [",
                           round(ci_death_diff_ss[1], 0),", ",
                           round(ci_death_diff_ss[2], 0), "]"),
                    paste0(round(ci_death_diff_pg[3], 0), " [",
                           round(ci_death_diff_pg[1], 0),", ",
                           round(ci_death_diff_pg[2], 0), "]"),
                    paste0(round(ci_death_diff_all[3], 0), " [",
                           round(ci_death_diff_all[1], 0),", ",
                           round(ci_death_diff_all[2], 0), "]"))

mean_oddeath_diff_nalox <- mean(m_outcomes_psa[, "oddeath_diff_nalox"])
ci_oddeath_diff_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_ss <- mean(m_outcomes_psa[, "oddeath_diff_ss"])
ci_oddeath_diff_ss <- quantile(m_outcomes_psa[, "oddeath_diff_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_pg <- mean(m_outcomes_psa[, "oddeath_diff_pg"])
ci_oddeath_diff_pg <- quantile(m_outcomes_psa[, "oddeath_diff_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_all <- mean(m_outcomes_psa[, "oddeath_diff_all"])
ci_oddeath_diff_all <- quantile(m_outcomes_psa[, "oddeath_diff_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_oddeath_diff <- c(paste0(round(ci_oddeath_diff_nalox[3], 0), " [",
                           round(ci_oddeath_diff_nalox[1], 0),", ",
                           round(ci_oddeath_diff_nalox[2], 0), "]"),
                    paste0(round(ci_oddeath_diff_ss[3], 0), " [",
                           round(ci_oddeath_diff_ss[1], 0),", ",
                           round(ci_oddeath_diff_ss[2], 0), "]"),
                    paste0(round(ci_oddeath_diff_pg[3], 0), " [",
                           round(ci_oddeath_diff_pg[1], 0),", ",
                           round(ci_oddeath_diff_pg[2], 0), "]"),
                    paste0(round(ci_oddeath_diff_all[3], 0), " [",
                           round(ci_oddeath_diff_all[1], 0),", ",
                           round(ci_oddeath_diff_all[2], 0), "]"))


effects_tbl <- cbind.data.frame(mod_names[-1], mean_cost_diff,
                                mean_death_diff, mean_oddeath_diff)

effects_tbl |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         mean_cost_diff = "Discounted Net Present \n Costs in Millions ",
                         mean_death_diff = "Deaths ",
                         mean_oddeath_diff = "Opioid-related Overdose Deaths")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Change Compared to No Intervention
Median [95% Credible Intervals]

Interventions

Discounted Net Present
Costs in Millions

Deaths

Opioid-related Overdose Deaths

Naloxone

209 [107, 362]

625 [414, 926]

-10089 [-15009, -6697]

Safer Supply

4216 [3651, 4821]

-14180 [-17532, -11074]

-3221 [-4105, -2529]

Prescription Guidelines

-25455 [-31480, -19866]

14417 [8610, 20565]

-5458 [-11045, -1624]

All Interventions

-21058 [-27059, -15403]

960 [-5687, 7753]

-18544 [-28496, -11309]

Code
mean_cost_diff_per_nalox <- mean(m_outcomes_psa[, "cost_diff_per_nalox"])
ci_cost_diff_per_nalox <- quantile(m_outcomes_psa[, "cost_diff_per_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_per_ss <- mean(m_outcomes_psa[, "cost_diff_per_ss"])
ci_cost_diff_per_ss <- quantile(m_outcomes_psa[, "cost_diff_per_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_per_pg <- mean(m_outcomes_psa[, "cost_diff_per_pg"])
ci_cost_diff_per_pg <- quantile(m_outcomes_psa[, "cost_diff_per_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_per_all <- mean(m_outcomes_psa[, "cost_diff_per_all"])
ci_cost_diff_per_all <- quantile(m_outcomes_psa[, "cost_diff_per_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_cost_diff_per <- c(paste0(round(ci_cost_diff_per_nalox[3], 2), " [",
                           round(ci_cost_diff_per_nalox[1], 2),", ",
                           round(ci_cost_diff_per_nalox[2], 2), "]%"),
                    paste0(round(ci_cost_diff_per_ss[3], 2), " [",
                           round(ci_cost_diff_per_ss[1], 2),", ",
                           round(ci_cost_diff_per_ss[2], 2), "]%"),
                    paste0(round(ci_cost_diff_per_pg[3], 2), " [",
                           round(ci_cost_diff_per_pg[1], 2),", ",
                           round(ci_cost_diff_per_pg[2], 2), "]%"),
                    paste0(round(ci_cost_diff_per_all[3], 2), " [",
                           round(ci_cost_diff_per_all[1], 2),", ",
                           round(ci_cost_diff_per_all[2], 2), "]%"))


mean_death_diff_per_nalox <- mean(m_outcomes_psa[, "death_diff_per_nalox"])
ci_death_diff_per_nalox <- quantile(m_outcomes_psa[, "death_diff_per_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_per_ss <- mean(m_outcomes_psa[, "death_diff_per_ss"])
ci_death_diff_per_ss <- quantile(m_outcomes_psa[, "death_diff_per_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_per_pg <- mean(m_outcomes_psa[, "death_diff_per_pg"])
ci_death_diff_per_pg <- quantile(m_outcomes_psa[, "death_diff_per_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_per_all <- mean(m_outcomes_psa[, "death_diff_per_all"])
ci_death_diff_per_all <- quantile(m_outcomes_psa[, "death_diff_per_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_death_diff_per <- c(paste0(round(ci_death_diff_per_nalox[3], 2), " [",
                           round(ci_death_diff_per_nalox[1], 2),", ",
                           round(ci_death_diff_per_nalox[2], 2), "]%"),
                    paste0(round(ci_death_diff_per_ss[3], 2), " [",
                           round(ci_death_diff_per_ss[1], 2),", ",
                           round(ci_death_diff_per_ss[2], 2), "]%"),
                    paste0(round(ci_death_diff_per_pg[3], 2), " [",
                           round(ci_death_diff_per_pg[1], 2),", ",
                           round(ci_death_diff_per_pg[2], 2), "]%"),
                    paste0(round(ci_death_diff_per_all[3], 2), " [",
                           round(ci_death_diff_per_all[1], 2),", ",
                           round(ci_death_diff_per_all[2], 2), "]%"))

mean_oddeath_diff_per_nalox <- mean(m_outcomes_psa[, "oddeath_diff_per_nalox"])
ci_oddeath_diff_per_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_per_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_per_ss <- mean(m_outcomes_psa[, "oddeath_diff_per_ss"])
ci_oddeath_diff_per_ss <- quantile(m_outcomes_psa[, "oddeath_diff_per_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_per_pg <- mean(m_outcomes_psa[, "oddeath_diff_per_pg"])
ci_oddeath_diff_per_pg <- quantile(m_outcomes_psa[, "oddeath_diff_per_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_per_all <- mean(m_outcomes_psa[, "oddeath_diff_per_all"])
ci_oddeath_diff_per_all <- quantile(m_outcomes_psa[, "oddeath_diff_per_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_oddeath_diff_per <- c(paste0(round(ci_oddeath_diff_per_nalox[3], 2), " [",
                           round(ci_oddeath_diff_per_nalox[1], 2),", ",
                           round(ci_oddeath_diff_per_nalox[2], 2), "]%"),
                    paste0(round(ci_oddeath_diff_per_ss[3], 2), " [",
                           round(ci_oddeath_diff_per_ss[1], 2),", ",
                           round(ci_oddeath_diff_per_ss[2], 2), "]%"),
                    paste0(round(ci_oddeath_diff_per_pg[3], 2), " [",
                           round(ci_oddeath_diff_per_pg[1], 2),", ",
                           round(ci_oddeath_diff_per_pg[2], 2), "]%"),
                    paste0(round(ci_oddeath_diff_per_all[3], 2), " [",
                           round(ci_oddeath_diff_per_all[1], 2),", ",
                           round(ci_oddeath_diff_per_all[2], 2), "]%"))


effects_tbl_per <- cbind.data.frame(mod_names[-1], mean_cost_diff_per,
                                mean_death_diff_per, mean_oddeath_diff_per)

effects_tbl_per |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Percent Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         mean_cost_diff_per = "Discounted Net Present Costs",
                         mean_death_diff_per = "Deaths ",
                         mean_oddeath_diff_per = "Opioid-related Overdose Deaths")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Percent Change Compared to No Intervention
Median [95% Credible Intervals]

Interventions

Discounted Net Present Costs

Deaths

Opioid-related Overdose Deaths

Naloxone

0.03 [0.01, 0.05]%

0.01 [0.01, 0.02]%

-3.46 [-4.13, -3.13]%

Safer Supply

0.54 [0.46, 0.61]%

-0.31 [-0.38, -0.24]%

-1.11 [-1.97, -0.66]%

Prescription Guidelines

-3.23 [-3.94, -2.55]%

0.32 [0.19, 0.45]%

-1.81 [-2.76, -0.93]%

All Interventions

-2.67 [-3.39, -1.98]%

0.02 [-0.13, 0.17]%

-6.43 [-7.16, -5.58]%

Code
tbl_effect_plot1 <- cbind.data.frame(value = c("lb", "ub", "median"),
                               ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
                               ci_cost_diff_per_pg, ci_cost_diff_per_all,
                               ci_death_diff_per_nalox, ci_death_diff_per_ss,
                               ci_death_diff_per_pg, ci_death_diff_per_all,
                               ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
                               ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>% 
  pivot_longer(2:13, names_to = "grp", values_to = "per_diff") %>% 
  pivot_wider(names_from = "value", values_from = "per_diff") %>% 
  separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>% 
  select(-c(type1, type2, type3)) %>% 
  mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
                               ifelse(Intervention == "ss", "Safer Supply",
                                      ifelse(Intervention == "pg", "Prescription Guidelines", 
                                             ifelse(Intervention == "all", "All Interventions", Intervention)))))
  
tbl_effect_plot1$Intervention <- factor(tbl_effect_plot1$Intervention, levels = mod_names, labels = mod_names)

ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = median))+
  geom_jitter(aes(color = Intervention, shape = Intervention), size = 2, position = position_dodge(0.25)) +
    # geom_errorbar(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), width = 0.3, position = position_dodge(0.1)) +
    geom_pointrange(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), position = position_dodge(0.25)) + 
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16, 17, 18)) +
  # scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","Opioid-related Overdose Deaths"))+
    theme_minimal())
Code
tbl_effect_plot2 <- cbind.data.frame(value = c("cost_lb", "cost_ub", "cost_median"),
                                     ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
                                     ci_cost_diff_per_pg, ci_cost_diff_per_all) %>% 
  pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
   pivot_wider(names_from = "value", values_from = "per_diff") %>% 
  separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>% 
  select(-c(type1, type2, type3, group)) %>% 
  full_join(., cbind.data.frame(value = c("oddeath_lb", "oddeath_ub", "oddeath_median"),
                                     ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
                                     ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>% 
              pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
              pivot_wider(names_from = "value", values_from = "per_diff") %>% 
              separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>% 
              select(-c(type1, type2, type3, group)), by = "Intervention") %>% 
  mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
                               ifelse(Intervention == "ss", "Safer Supply",
                                      ifelse(Intervention == "pg", "Prescription Guidelines", 
                                             ifelse(Intervention == "all", "All Interventions", Intervention)))))
  
  tbl_effect_plot2$Intervention <- factor(tbl_effect_plot2$Intervention, levels = mod_names, labels = mod_names)


ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_median, y = oddeath_median)) +
  geom_point(aes(color = Intervention), size = 1) +
    geom_pointrange(aes(ymin = oddeath_lb, ymax = oddeath_ub, color = Intervention)) +
  geom_errorbarh(aes(xmax = cost_lb, xmin = cost_ub,color = Intervention),  height = 0) +
  scale_color_brewer(palette = "Dark2") +
  # scale_shape_manual(values = c(15, 16, 17, 18)) +
    geom_hline(yintercept = 0, linewidth = 0.25) +
  geom_vline(xintercept = 0, linewidth = 0.25) +
  xlab("Change in Costs Compared to No Intervention (%)") +
  ylab("Change in Opioid-related Overdose Deaths \n Compared to No Intervention (%)") +
  theme_minimal())

The results show the 3 interventions individually and all together reduced the total number of opioid-related deaths compared to no interventions, with all interventions combined resulting in -6.413539 [-5.58, -7.16] % difference compared to no intervention scenario. All interventions and prescription guidelines resulted in reduction costs as well, however Naloxone and Safer Supply increased costs slightly. Implementing all interventions reduces both costs and opioid-related OD deaths.
Regarding overall deaths though, except for safer supply, all other scenarios either increase overall deaths or result in inconclusive results.

3.3 Results for the secondary outcomes (point estimates)

Code
v_yrs_prev <- c(2015:2029)
yrs_prev <- length(v_yrs_prev)

prev_fun <- function(mod, i, v_states){
  
  prop_uncalib <- data.frame(matrix(NA,byrow = T, nrow = 12, 
                                             ncol = yrs_prev,
                                             list(NULL,
                                                  paste0("prop_",
                                                         str_sub(v_yrs_prev, 3, 4)))))
  
  tot_pop_uncalib_tbl <- data.frame(matrix(NA,byrow = T, nrow = 12, 
                                             ncol = yrs_prev,
                                 list(NULL,
                                      paste0("tot_pop_",
                                             str_sub(v_yrs_prev, 3, 4)))))
  
  for (j in 1:yrs_prev) {
    yr <- (v_yrs_prev[1] - 1) + j
    
    
    if (length(v_states) > 1)
      {
      prop_uncalib[, j] <- assign(
        paste0("prop_", str_sub(yr, 3, 4)),
        rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             v_states]) /
          rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
        )
    } else {
      prop_uncalib[, j] <- assign(
      paste0("prop_", str_sub(yr, 3, 4)),
      mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             v_states]) /
      rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
    }
    
    
  
  tot_pop_uncalib_tbl[, j] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  
  }
  
  # final table for monthly prevalence over 15 years
  prop_uncalib_tbl <- prop_uncalib %>% 
    pivot_longer(1:15, names_to = "grp", values_to = "value") %>% 
    arrange(grp) %>% 
    mutate(year = rep(v_yrs_prev,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_uncalib_tbl %>% 
              pivot_longer(1:15, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop)) %>% 
    mutate(intervention = mod_names[i])
  
  # table with weighted yearly prevalence
  prop_uncalib_wei_mean <- prop_uncalib_tbl %>% 
    group_by(year) %>% 
    summarise(value = weighted.mean(value, tot_pop_val)) %>% 
    ungroup() %>% 
    mutate(intervention = mod_names[i]) 
  
  return (list(prop_uncalib_tbl = prop_uncalib_tbl,
               prop_uncalib_wei_mean = prop_uncalib_wei_mean))
}

# Prevalence of severe brain injury

noint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
nalx_prev_sevbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
ss_prev_sevbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
pg_prev_sevbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
allint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean


prev_sevbi_yr_tbl <- noint_prev_sevbi_yr_tbl %>% 
  bind_rows(., nalx_prev_sevbi_yr_tbl, ss_prev_sevbi_yr_tbl,
            pg_prev_sevbi_yr_tbl, allint_prev_sevbi_yr_tbl)

Weighted Yearly Prevalence of Severe Brain Injury

Code
ggplotly(ggplot() +
  geom_point(data = prev_sevbi_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
    theme_minimal())
Code
# Prevalence of moderate Brain injury

noint_prev_modbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
nalx_prev_modbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
ss_prev_modbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
pg_prev_modbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
allint_prev_modbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean

prev_modbi_yr_tbl <- noint_prev_modbi_yr_tbl %>% 
  bind_rows(., nalx_prev_modbi_yr_tbl, ss_prev_modbi_yr_tbl,
            pg_prev_modbi_yr_tbl, allint_prev_modbi_yr_tbl)

Weighted Yearly Prevalence of Moderate Brain Injury

Code
ggplotly(ggplot() +
  geom_point(data = prev_modbi_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())
Code
# Prevalence of substance use treatment

v_states_prev_tx <- c("BS_DETOX", "BS_OAT_INI", "BS_OAT_MAINT",
                          "BS_OAT_SS", "BS_SS", "BR_OAT_INI", "BR_OAT_MAINT",
                          "BR_OAT_SS", "BR_SS")
noint_prev_tx_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
nalx_prev_tx_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
ss_prev_tx_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
pg_prev_tx_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
allint_prev_tx_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean

prev_tx_yr_tbl <- noint_prev_tx_yr_tbl %>% 
  bind_rows(., nalx_prev_tx_yr_tbl, ss_prev_tx_yr_tbl,
            pg_prev_tx_yr_tbl, allint_prev_tx_yr_tbl)

Weighted Yearly Prevalence of substance use treatment

Code
ggplotly(ggplot() +
  geom_point(data = prev_tx_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())
Code
# Prevalence of prescription opioid use

prop_opioids_rx_target_tbl <- calib_target_tbl %>% 
              filter(group == "prev_on_opioidsrx") %>% 
              select(year, target) %>% 
              rename(value = target) %>% 
              mutate(year_mon = paste(year, month.abb[7], sep = "_"))
  
v_states_prev_rx_opd <- c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")

noint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
nalx_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
ss_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
pg_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
allint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
  

prev_rx_opd_yr_tbl <- noint_prev_rx_opd_yr_tbl %>% 
  bind_rows(., nalx_prev_rx_opd_yr_tbl, ss_prev_rx_opd_yr_tbl,
            pg_prev_rx_opd_yr_tbl, allint_prev_rx_opd_yr_tbl) %>% 
    bind_rows(., prop_opioids_rx_target_tbl %>% 
              select(year, value) %>% 
              mutate(intervention = "Target"))

Weighted Yearly Prevalence of prescription opioid use

Code
ggplotly(ggplot() +
  geom_point(data = prev_rx_opd_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())
Code
# Prevalence of illicit use

v_states_illicit <- c("BI_ILLICIT", "BR_ILLICIT")
noint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
nalx_prev_illicituse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
ss_prev_illicituse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
pg_prev_illicituse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
allint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean

prev_illicituse_yr_tbl <- noint_prev_illicituse_yr_tbl %>% 
  bind_rows(., nalx_prev_illicituse_yr_tbl, ss_prev_illicituse_yr_tbl,
            pg_prev_illicituse_yr_tbl, allint_prev_illicituse_yr_tbl)

Weighted Yearly Prevalence of illicit opioid use

Code
ggplotly(ggplot() +
  geom_point(data = prev_illicituse_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())
Code
# Prevalence of misuse

noint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
nalx_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
ss_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
pg_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
allint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean

prev_rxmisuse_yr_tbl <- noint_prev_rxmisuse_yr_tbl %>% 
  bind_rows(., nalx_prev_rxmisuse_yr_tbl, ss_prev_rxmisuse_yr_tbl,
            pg_prev_rxmisuse_yr_tbl, allint_prev_rxmisuse_yr_tbl)

Weighted Yearly Prevalence of rx opioid misuse

Code
ggplotly(ggplot() +
  geom_point(data = prev_rxmisuse_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())
Code
# Prevalence of overdose
v_states_od <- c("BO_OD_RX", "BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_od_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
nalx_prev_od_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
ss_prev_od_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
pg_prev_od_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
allint_prev_od_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_od)$prop_uncalib_wei_mean

prev_od_yr_tbl <- noint_prev_od_yr_tbl %>% 
  bind_rows(., nalx_prev_od_yr_tbl, ss_prev_od_yr_tbl,
            pg_prev_od_yr_tbl, allint_prev_od_yr_tbl)

Weighted Yearly Prevalence of opioid-related overdose

Code
ggplotly(ggplot() +
  geom_point(data = prev_od_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + theme_minimal())

The results show that none of the interventions affected prevalence of brain injury except for Safer Supply. Safer supply also resulted in an increase of prevalence in substance use treatment (which is expected since it is part of substance use treatment). Prescription Guidelines and All interventions resulted in a reduction in prevalence of prescription opioid use and prescription opioid misuse. No notable differences were observed for prevalence of opioid illicit drug use and opioid-related overdose.

4 Discussion and next steps

What is the relationship between your findings and the wider literature or policy landscape?

This analysis shows the promising effects of these interventions in reducing the direct harms of opioid crisis by reducing opioid-related OD deaths. From a healthcare perspective, when all interventions are implemented simultaneously they would reduce costs while saving lives. However, the results regarding overall deaths are not as promising, we are interested in examining the pathways that led to that increase in overall deaths.
One of the strengths of this analysis is that the model distinguished between chronic noncancer and cancer pain and included palliative care. It also included brain injuries as one of the often neglected sequalea of overdose. The uncalibrated model estimated values that were similar to target values (a strength of the model), however, it seemed that the calibration did not make a significant difference (a limitation). The calibration might not have worked due to non-identifiability problem since we have many parameters that we were trying to calibrate and few calibration targets. However that is due to the lack of data availability. This has been the main limitation of this project, there is a lack of available for the different parameters, many of which were based on scientific assumptions and knowledge of experts in the field. Therefore, as a next step, we are interested in 1) repeating the analysis for only Ontario since they have more data available so we can better calibrate the model parameters, 2) conduct a value of information analysis to identify the parameters that we need to collect data on, 3) use British Colombia data to try to further understand the pathways that led to increase in overall deaths and repeat the analysis for this province similar to Onatrio and 4) conduct the same analysis with societal costs in addition to healthcare costs.

5 References

6 Appendix 1: Tables of the parameters used in the analysis

6.1 Initial Population counts in each state

Code
flextable(ini_states_tbl_org %>% 
            select(`Variable `, basevalue)) |>
  autofit()

Variable 

basevalue

BN_PN

12,329,264

BN_ACUTE

77,121

BN_CHRONIC

6,000,000

BN_CANCER

163,129

BN_OTHER

6,003,590

BN_PALLIATIVE

4,003

BPO_ACUTE

216,300

BPO_CHRONIC

3,530,940

BPO_CANCER

123,666

BPO_OTHER

594,870

BPO_PALLIATIVE

14,158

BPO_MISUSE

274,426

BI_ILLICIT

290,063

BS_DETOX

7,500

BS_OAT_INI

2,461

BS_OAT_MAINT

72,573

BS_OAT_SS

0

BS_SS

0

BO_OD_ILLICIT

0

BO_OD_RX

0

BO_MOD_BI

0

BO_SEVERE_BI

0

BO_SEVERE_BI_OUT

0

BO_OD_DEATH

0

BO_DEATH

0

BR_ILLICIT

0

BR_OAT_INI

0

BR_OAT_MAINT

0

BR_OAT_SS

0

BR_SS

0

BR_OD_ILLICIT

0

6.2 Population increase every month

Code
flextable(pop_increase_tbl %>% 
            mutate(Year = as.character(Year)) %>% 
            rename("Increase per month" = increase_per_month)) |>
  autofit()

Year

Increase per month

2015

18,637

2016

27,882

2017

32,458

2018

38,812

2019

41,247

2020

31,231

2021

19,700

2022

54,525

2023

48,524

2024

51,161

2025

50,779

2026

50,018

2027

49,005

2028

48,006

2029

47,094

6.3 Transition probabilities

Code
flextable(trans_prob_tbl_new %>%
            filter(basevalue != 0) %>% 
            rename("From - state" = var_from,
                   "To - state" = var_to,
                   Basevalue = basevalue)) |>
  autofit()

From - state

To - state

Basevalue

BN_PN

BN_PN

999.00000000

BN_PN

BN_ACUTE

0.00165257

BN_PN

BN_CANCER

0.00050000

BN_PN

BN_OTHER

0.00069143

BN_PN

BN_PALLIATIVE

0.00030770

BN_PN

BPO_ACUTE

0.00110106

BN_PN

BPO_MISUSE

0.00220011

BN_PN

BI_ILLICIT

0.00002940

BN_PN

BO_DEATH

0.00032167

BN_ACUTE

BN_PN

999.00000000

BN_ACUTE

BN_CHRONIC

0.60000000

BN_ACUTE

BPO_MISUSE

0.00004680

BN_ACUTE

BI_ILLICIT

0.00001040

BN_ACUTE

BO_DEATH

0.00053201

BN_CHRONIC

BN_PN

0.00010000

BN_CHRONIC

BN_CHRONIC

999.00000000

BN_CHRONIC

BN_CANCER

0.00042924

BN_CHRONIC

BN_PALLIATIVE

0.00030770

BN_CHRONIC

BPO_CHRONIC

0.01000000

BN_CHRONIC

BPO_MISUSE

0.00728000

BN_CHRONIC

BI_ILLICIT

0.00010000

BN_CHRONIC

BO_DEATH

0.00040000

BN_CANCER

BN_PN

0.00010000

BN_CANCER

BN_CHRONIC

0.00003130

BN_CANCER

BN_CANCER

999.00000000

BN_CANCER

BN_PALLIATIVE

0.00066950

BN_CANCER

BPO_CANCER

0.55000000

BN_CANCER

BO_DEATH

0.01390000

BN_OTHER

BN_PN

0.00053731

BN_OTHER

BN_OTHER

999.00000000

BN_OTHER

BPO_OTHER

0.00349866

BN_OTHER

BO_DEATH

0.00060864

BN_PALLIATIVE

BN_PALLIATIVE

0.18500000

BN_PALLIATIVE

BPO_PALLIATIVE

0.44000000

BN_PALLIATIVE

BO_DEATH

999.00000000

BPO_ACUTE

BN_PN

999.00000000

BPO_ACUTE

BN_CHRONIC

0.60000000

BPO_ACUTE

BPO_CHRONIC

0.06200000

BPO_ACUTE

BPO_MISUSE

0.00600000

BPO_ACUTE

BO_OD_RX

0.00743025

BPO_ACUTE

BO_DEATH

0.00053201

BPO_CHRONIC

BN_CHRONIC

0.01750000

BPO_CHRONIC

BPO_CHRONIC

999.00000000

BPO_CHRONIC

BPO_MISUSE

0.00470310

BPO_CHRONIC

BO_OD_RX

0.00160398

BPO_CHRONIC

BO_DEATH

0.00040000

BPO_CANCER

BN_CANCER

0.00874161

BPO_CANCER

BPO_CHRONIC

0.00177217

BPO_CANCER

BPO_CANCER

999.00000000

BPO_CANCER

BPO_PALLIATIVE

0.00100425

BPO_CANCER

BPO_MISUSE

0.00470310

BPO_CANCER

BO_OD_RX

0.00097898

BPO_CANCER

BO_DEATH

0.01390000

BPO_OTHER

BN_OTHER

0.02825010

BPO_OTHER

BPO_OTHER

999.00000000

BPO_OTHER

BPO_MISUSE

0.00470310

BPO_OTHER

BO_OD_RX

0.00097898

BPO_OTHER

BO_DEATH

0.00040000

BPO_PALLIATIVE

BPO_PALLIATIVE

999.00000000

BPO_PALLIATIVE

BO_DEATH

0.38250000

BPO_MISUSE

BN_PN

0.01000000

BPO_MISUSE

BPO_MISUSE

999.00000000

BPO_MISUSE

BI_ILLICIT

0.00426953

BPO_MISUSE

BS_DETOX

0.00100000

BPO_MISUSE

BS_OAT_INI

0.00100000

BPO_MISUSE

BO_OD_RX

0.00097898

BPO_MISUSE

BO_DEATH

0.00064333

BI_ILLICIT

BN_PN

0.00001000

BI_ILLICIT

BN_CANCER

0.00064386

BI_ILLICIT

BI_ILLICIT

999.00000000

BI_ILLICIT

BS_DETOX

0.01000000

BI_ILLICIT

BS_OAT_INI

0.00700000

BI_ILLICIT

BO_OD_ILLICIT

0.01320000

BI_ILLICIT

BO_DEATH

0.00080801

BS_DETOX

BN_PN

999.00000000

BS_DETOX

BI_ILLICIT

0.22181118

BS_DETOX

BS_OAT_INI

0.82000000

BS_DETOX

BO_OD_ILLICIT

0.01431324

BS_DETOX

BO_DEATH

0.00309346

BS_OAT_INI

BI_ILLICIT

999.00000000

BS_OAT_INI

BS_OAT_MAINT

0.91700000

BS_OAT_INI

BO_OD_ILLICIT

0.04050000

BS_OAT_INI

BO_DEATH

0.00150000

BS_OAT_MAINT

BN_PN

0.00460000

BS_OAT_MAINT

BI_ILLICIT

999.00000000

BS_OAT_MAINT

BS_OAT_MAINT

0.71000000

BS_OAT_MAINT

BO_OD_ILLICIT

0.00420000

BS_OAT_MAINT

BO_DEATH

0.00038326

BS_OAT_SS

BS_OAT_SS

999.00000000

BS_SS

BS_SS

999.00000000

BO_OD_ILLICIT

BI_ILLICIT

999.00000000

BO_OD_ILLICIT

BS_DETOX

0.26800000

BO_OD_ILLICIT

BS_OAT_INI

0.05871300

BO_OD_ILLICIT

BO_MOD_BI

0.03000000

BO_OD_ILLICIT

BO_SEVERE_BI

0.00086172

BO_OD_ILLICIT

BO_OD_DEATH

0.01900000

BO_OD_RX

BPO_MISUSE

999.00000000

BO_OD_RX

BS_DETOX

0.26800000

BO_OD_RX

BS_OAT_INI

0.05871300

BO_OD_RX

BO_MOD_BI

0.03000000

BO_OD_RX

BO_SEVERE_BI

0.00086172

BO_OD_RX

BO_OD_DEATH

0.01380000

BO_MOD_BI

BO_DEATH

0.00874161

BO_MOD_BI

BR_ILLICIT

999.00000000

BO_MOD_BI

BR_OAT_INI

0.15100000

BO_SEVERE_BI

BO_SEVERE_BI_OUT

999.00000000

BO_SEVERE_BI

BO_DEATH

0.67000000

BO_SEVERE_BI_OUT

BO_SEVERE_BI_OUT

999.00000000

BO_SEVERE_BI_OUT

BO_DEATH

0.00090000

BO_OD_DEATH

BO_OD_DEATH

1.00000000

BO_DEATH

BO_DEATH

1.00000000

BR_ILLICIT

BO_DEATH

0.01740000

BR_ILLICIT

BR_ILLICIT

999.00000000

BR_ILLICIT

BR_OAT_INI

0.02712532

BR_ILLICIT

BR_OD_ILLICIT

0.03650000

BR_OAT_INI

BO_DEATH

0.01740000

BR_OAT_INI

BR_ILLICIT

999.00000000

BR_OAT_INI

BR_OAT_MAINT

0.84075000

BR_OAT_INI

BR_OD_ILLICIT

0.01285868

BR_OAT_MAINT

BO_DEATH

0.00109940

BR_OAT_MAINT

BR_ILLICIT

999.00000000

BR_OAT_MAINT

BR_OAT_MAINT

0.84075000

BR_OAT_MAINT

BR_OD_ILLICIT

0.01820000

BR_OAT_SS

BR_OAT_SS

999.00000000

BR_SS

BR_SS

999.00000000

BR_OD_ILLICIT

BO_MOD_BI

999.00000000

BR_OD_ILLICIT

BO_SEVERE_BI

0.03000000

BR_OD_ILLICIT

BO_OD_DEATH

0.07700000

6.4 Healthcare costs for all states used in the analysis

Code
flextable(costs_tbl) |>
  autofit()

Variable 

Cost

BN_PN

0.00

BN_ACUTE

200.34

BN_CHRONIC

70.93

BN_CANCER

580.18

BN_OTHER

70.93

BN_PALLIATIVE

5,405.14

BPO_ACUTE

15,002.06

BPO_CHRONIC

103.05

BPO_CANCER

2,186.73

BPO_OTHER

103.05

BPO_PALLIATIVE

2,904.79

BPO_MISUSE

912.27

BI_ILLICIT

97.10

BS_DETOX

4,550.00

BS_OAT_INI

1,500.88

BS_OAT_MAINT

1,393.88

BS_OAT_SS

1,432.77

BS_SS

1,780.69

BO_OD_RX

974.48

BO_OD_ILLICIT

974.48

BO_MOD_BI

6,825.46

BO_SEVERE_BI

46,749.90

BO_SEVERE_BI_OUT

1,295.96

BO_OD_DEATH

0.00

BO_DEATH

0.00

BR_ILLICIT

97.10

BR_OAT_INI

1,500.88

BR_OAT_MAINT

1,393.88

BR_OAT_SS

1,432.77

BR_SS

1,780.69

BR_OD_ILLICIT

97.10

7 Appendix 2: Comparison between Uncalibrated model and Calibration Targets

Code
library(here)
library(RColorBrewer)
source(here("02_scripts/01_fun_data.R"))
calib_target_tbl <- read_excel(here("01_data/markov_model_parameters_preprior.xlsx"),
                               sheet = "calibration_targets")
theme_set(theme_minimal())

7.1 Prevalence of Rx opioids use

Code
# prevalence of people on prescribed opioids

v_yrs_prev_rx <- c(2015:2018)
yrs_prev_rx <- length(v_yrs_prev_rx)

prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
                                             nrow = 12,
                                             ncol = yrs_prev_rx,
                                             list(NULL,
                                                  paste0("prop_opioids_rx_",
                                                         str_sub(v_yrs_prev_rx, 3, 4)))))
tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
                                         nrow = 12,
                                         ncol = yrs_prev_rx,
                                 list(NULL,
                                      paste0("tot_pop_",
                                             str_sub(v_yrs_prev_rx, 3, 4)))))

for (i in 1:yrs_prev_rx) {
  yr <- v_yrs_prev_rx[1] + i
  prop_opioids_rx_uncalib[, i] <- assign(
    paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")]) /
      rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
  )
  
  tot_pop_uncalib_tbl[, i] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  }

prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>% 
  pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>% 
  arrange(grp) %>% 
  mutate(year = rep(v_yrs_prev_rx,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_uncalib_tbl %>% 
              pivot_longer(1:4, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop))

prop_opioids_rx_target_tbl <- calib_target_tbl %>% 
              filter(group == "prev_on_opioidsrx") %>% 
              select(year, target, group) %>% 
              rename(grp = group,
                     target_val = target) %>% 
              mutate(year_mon = paste(year, month.abb[7], sep = "_"))

prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>% 
  group_by(year) %>% 
  summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
  ungroup() %>% 
  mutate(grp = "Model") %>% 
  bind_rows(., prop_opioids_rx_target_tbl %>% 
              select(year, target_val) %>% 
              mutate(grp = "Target"))

prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
                                              levels = paste(rep(v_yrs_prev_rx,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(v_yrs_prev_rx,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
                                              levels = paste(rep(v_yrs_prev_rx,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(v_yrs_prev_rx,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))

p1 <- ggplot() +
  geom_point(data = prop_opioids_rx_uncalib_tbl,
             aes(x = year_mon, y = target_val,
                 color = factor(year), fill = factor(year))) +
  geom_point(data = prop_opioids_rx_target_tbl,
             aes(x = year_mon, y = target_val),
             color = "red") +
  geom_point(data = prop_opioids_rx_uncalib_tbl %>% 
               group_by(year) %>% 
               summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(year_mon = prop_opioids_rx_target_tbl$year_mon),
             aes(x = year_mon, y = target_val),
             color = "blue")+
  xlab("Year_Month") + ylab("Prevalence of prescription opioid use") +
  scale_fill_discrete(name = "year")+
  scale_color_discrete(name = "year")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

plotly::ggplotly(p1)

The red points represent the observed target proportions of Rx opioids from CIHI report, the blue points represent an annual average (weighted by population during that month), the other coloured points represent monthly prevalence predicted from the model

Code
p2 <- ggplot() +
  geom_point(data = prop_opioids_rx_uncalib_wei_mean,
             aes(x = year, y = target_val, color = factor(grp),
                 fill = factor(grp))) +
  xlab("Year") + ylab("Prevalence of prescription opioid use")+
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
  # ylim(0,0.2) +
  # ggtitle("Comparison between uncalibrated annual weighted average and
  #         observed targets of prevalence of Rx opioids use")

plotly::ggplotly(p2)

7.2 Deaths and OD-Deaths

Code
################ Deaths 
# number of deaths
v_yr1_deaths <- c(2016:2020)
yrs_deaths <- length(v_yr1_deaths)

num_deaths_uncalib <- rep(NA, yrs_deaths)
for (i in 1:yrs_deaths){
  yr <- (v_yr1_deaths[1] - 1) + i
  num_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                       year_mon_cycle_tbl$mon == 12] + 1,
                                            "BO_DEATH"] -
    mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_DEATH"]
  }
  
# number of OD-deaths
v_yr1_oddeaths <- c(2016:2021)
yrs_oddeaths <- length(v_yr1_oddeaths)

num_od_deaths_uncalib <- rep(NA, yrs_oddeaths)
for (i in 1:yrs_oddeaths){
  yr <- (v_yr1_oddeaths[1] - 1) + i
  num_od_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
}


deaths_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% c("total_deaths", "total_od_deaths")) %>% 
  mutate(target = ifelse(target == "total_deaths", "Total deaths",
                         ifelse(target == "total_od_deaths",
                                "Total opioid-related overdose deaths", NA)),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_deaths_uncalib) %>% 
              mutate(target = "Total deaths") %>% 
              bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
                          mutate(target = "Total opioid-related overdose deaths")) %>% 
              mutate(year = c(v_yr1_deaths, v_yr1_oddeaths),
                     group = "Model"))

p3 <- ggplot() +
  geom_point(data = deaths_target_uncalib_tbl,
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Deaths") +
  scale_color_brewer(palette = "Dark2", name = "")  +
  scale_fill_brewer(palette = "Dark2", name = "") +
  facet_wrap(~target, scales = "free") +
  theme(legend.position="bottom")

plotly::ggplotly(p3)
Code
plotly::ggplotly(ggplot() +
  geom_point(data = deaths_target_uncalib_tbl %>% 
               filter(target == "Total deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Deaths") +
  scale_color_brewer(palette = "Dark2", name = "")  +
  scale_fill_brewer(palette = "Dark2", name = ""))
Code
  # facet_wrap(~target, scales = "free") + theme(legend.position="bottom")

plotly::ggplotly(ggplot() +
  geom_point(data = deaths_target_uncalib_tbl %>%
               filter(target == "Total opioid-related overdose deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total opioid-related overdose deaths") +
  scale_color_brewer(palette = "Dark2", name = "")  +
  scale_fill_brewer(palette = "Dark2", name = ""))

7.3 OAT

Code
# number of oat
v_yr1_oat <- c(2018:2021)
yrs_oat <- length(v_yr1_oat)
  
num_oat_uncalib <- rep(NA, yrs_oat)

for (i in 1:yrs_oat){
  yr <- (v_yr1_oat - 1) + i 
    num_oat_uncalib[i] <- sum(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                        year_mon_cycle_tbl$mon == 6] + 1,
                             c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}


oat_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% "total_oat") %>% 
  mutate(target = ifelse(target == "total_oat",
                                "Total OAT", NA),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
                          mutate(target = "Total OAT",
                                 year = v_yr1_oat,
                                 group = "Model"))

p3 <- ggplot() +
  geom_point(data = oat_target_uncalib_tbl,
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total") +
  scale_color_brewer(palette = "Dark2", name = "")  +
  scale_fill_brewer(palette = "Dark2", name = "") +
  theme(legend.position="bottom")



plotly::ggplotly(p3)

7.4 Trace plot

Code
# build-in color palette
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
                                        c("#084C61", "#DB504A",
                                                   "#E3B505", "#4F6D7A", 
                                                      "#56A3A6"))
names(colors_pal) <- NULL


trace_tbl_inc <- as_tibble(mod_basecase$m_M) %>% 
  mutate(cycle_num = 0:n_cycles) %>% 
  pivot_longer(1:n_states, names_to = "state", values_to = "count")


trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
                      trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_OD_DEATH" &
                      trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_od_deaths

trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
                      trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_MOD_BI" &
                      trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_mod_bi

trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
                      trace_tbl_inc$cycle_num == 180] <- trace_tbl_inc$count[trace_tbl_inc$state == "BO_SEVERE_BI" &
                      trace_tbl_inc$cycle_num == 180] + mod_basecase$extra_sev_bi

trace_tbl_inc$state <- factor(trace_tbl_inc$state,
                              levels = v_state_names,
                              labels = v_state_names)
p4 <- ggplot(trace_tbl_inc, aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl_inc %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal)

plotly::ggplotly(p4)
Code
p5 <- ggplot(trace_tbl_inc %>% 
               filter(state == "BO_OD_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               mutate(year = c(2015:2029)),
             aes(x = year, y = count)) +
  geom_line() +
  ylab("Total opioid-related overdose deaths") + xlab("Year") 
# +
#   ggtitle("Total opioid-related overdose deaths over time")

plotly::ggplotly(p5)

8 Appendix 3: Comparison between MAP calibrated model, uncalibrated model and calibration Targets

8.1 Scenario 1

It takes into account fentanyl contamination

Code
library(here)

source(here("02_scripts/03a_fun_calib_td.R"))
theme_set(theme_minimal())

opt_params <- readRDS(here("01_data/map_calib_td.RDS"))

Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets

Code
mod_opt_map <- mod_calib(as.numeric(opt_params$par))

8.2 Prevalence of Rx opioids use

Code
# prevalence of people on prescribed opioids
v_yrlast_prev_rx <- c(2015:2018)
yrs_prev_rx <- length(v_yrlast_prev_rx)

prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
                                             nrow = 12,
                                             ncol = yrs_prev_rx,
                                             list(NULL,
                                                  paste0("prop_opioids_rx_",
                                                         str_sub(v_yrlast_prev_rx, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
                                         nrow = 12,
                                         ncol = yrs_prev_rx,
                                 list(NULL,
                                      paste0("tot_pop_",
                                             str_sub(v_yrlast_prev_rx, 3, 4)))))

for (i in 1:yrs_prev_rx) {
  yr <- (v_yrlast_prev_rx[1] - 1) + i
  prop_opioids_rx_uncalib[, i] <- assign(
    paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")]) /
      rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
  )
  
  tot_pop_uncalib_tbl[, i] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  
  prop_opioids_rx_calib[, i] <- assign(
    paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
    rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")]) /
      rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
  )
  
  tot_pop_calib_tbl[, i] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  }

prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>% 
  pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>% 
  arrange(grp) %>% 
  mutate(year = rep(2015:2018,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_uncalib_tbl %>% 
              pivot_longer(1:4, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop))


prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>% 
  pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>% 
  arrange(grp) %>% 
  mutate(year = rep(2015:2018,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_calib_tbl %>% 
              pivot_longer(1:4, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop))


prop_opioids_rx_target_tbl <- calib_target_tbl %>% 
              filter(group == "prev_on_opioidsrx") %>% 
              select(year, target, group) %>% 
              rename(grp = group,
                     target_val = target) %>% 
              mutate(year_mon = paste(year, month.abb[6], sep = "_"))

prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>% 
  group_by(year) %>% 
  summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
  ungroup() %>% 
  mutate(grp = "Uncalibrated Model") %>% 
  bind_rows(., prop_opioids_rx_target_tbl %>% 
              select(year, target_val) %>% 
              mutate(grp = "target")) %>% 
  bind_rows(., prop_opioids_rx_calib_tbl %>% 
              group_by(year) %>% 
              summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
              ungroup() %>% 
              mutate(grp = "Calibrated Model"))

prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))

prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>% 
  mutate(group = "Uncalibrated Model") %>% 
  bind_rows(.,prop_opioids_rx_calib_tbl %>% 
              mutate(group = "Calibrated Model"))

wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>% 
  mutate(grp = "Target") %>% 
  bind_rows(., prop_opioids_rx_uncalib_tbl %>% 
               group_by(year) %>% 
               summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
                      grp = "Uncalibrated Model")) %>% 
  bind_rows(., prop_opioids_rx_calib_tbl %>% 
               group_by(year) %>% 
               summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
                      grp = "Calibrated Model"))

prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
                                  labels = c("Calibrated Model", "Uncalibrated Model"),
                                  levels = c("Calibrated Model", "Uncalibrated Model"))

p1 <- ggplot() +
  geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
                                             color = group, fill = group)) +
  geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
             aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 1.5) +
  scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

plotly::ggplotly(p1)

8.2.1 Comparison between uncalibrated annual weighted average and observed targets of prevalence of Rx opioids use

Code
wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p2 <- ggplot() +
 geom_point(data = wprop_opioids_rx_tbl,
             aes(x = year_mon, y = target_val,
                 color = grp, fill = grp), size = 1.5) +
  xlab("Year") + ylab("Prevalence of Rx opioid use")+
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
  # ylim(0,0.2) +
plotly::ggplotly(p2)

8.3 Deaths and OD-Deaths

Code
################ Deaths 
# number of deaths

num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
  yr <- 2015 + i
  num_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                       year_mon_cycle_tbl$mon == 12] + 1,
                                            "BO_DEATH"] -
    mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_DEATH"]
  
  num_deaths_calib[i] <- mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                       year_mon_cycle_tbl$mon == 12] + 1,
                                            "BO_DEATH"] -
    mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_DEATH"]
  }
  
# number of OD-deaths
  
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
  yr <- 2015 + i 
  num_od_deaths_uncalib[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  num_od_deaths_calib[i] <- mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
}

deaths_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% c("total_deaths", "total_od_deaths")) %>% 
  mutate(target = ifelse(target == "total_deaths", "Total deaths",
                         ifelse(target == "total_od_deaths",
                                "Total OD-deaths", NA)),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_deaths_uncalib) %>% 
              mutate(target = "Total deaths") %>% 
              bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
                          mutate(target = "Total OD-deaths")) %>% 
              mutate(year = c(2016:2020, 2016:2021),
                     group = "Uncalibrated Model")) %>% 
   bind_rows(., data.frame(target_val = num_deaths_calib) %>% 
              mutate(target = "Total deaths") %>% 
              bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
                          mutate(target = "Total OD-deaths")) %>% 
              mutate(year = c(2016:2020, 2016:2021),
                     group = "Calibrated Model"))
  
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p3 <- ggplot() +
  geom_point(data = deaths_target_uncalib_tbl  %>% filter(target == "Total deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Deaths") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
# +
#   facet_wrap(~target, scales = "free") + theme(legend.position="bottom")



plotly::ggplotly(p3)
Code
p4 <- ggplot() +
  geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
# +
#   facet_wrap(~target, scales = "free") + theme(legend.position="bottom")



plotly::ggplotly(p4)

8.4 OAT

Code
# number of oat
  
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
  yr <- 2017 + i 
    num_oat_uncalib[i] <- sum(mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                        year_mon_cycle_tbl$mon == 6] + 1,
                             c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
    
    num_oat_calib[i] <- sum(mod_opt_map$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                        year_mon_cycle_tbl$mon == 6] + 1,
                             c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}


oat_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% c("total_oat")) %>% 
  mutate(target = ifelse(target == "total_oat",
                         "Total OAT", NA),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
                          mutate(target = "Total OAT",
                                 year = c(2018:2021),
                                 group = "Uncalibrated Model")) %>% 
   bind_rows(.,  data.frame(target_val = num_oat_calib) %>%
                          mutate(target = "Total OAT",
                                 year = c(2018:2021),
                                 group = "Calibrated Model"))

oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))


p5 <- ggplot() +
  geom_point(data = oat_target_uncalib_tbl,
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total OAT counts") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
# +
#   facet_wrap(~target, scales = "free") + theme(legend.position="bottom")



plotly::ggplotly(p5)
Code
cbind.data.frame(parameter = calib_param_tbl$var_params,
                 basevalue = calib_param_tbl$basevalue,
                 opt_map_nm = opt_params$par) %>% 
  filter(basevalue != opt_map_nm)
                      parameter  basevalue opt_map_nm
1          p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2      p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3   p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4       p__BS_DETOX__BI_ILLICIT 0.18181118 0.22181118
5       p__BS_DETOX__BS_OAT_INI 0.80000000 0.82000000
6 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.71000000

8.5 Scenario 2

It doesn’t take into account fentanyl contamination

Code
library(here)

source(here("02_scripts/03b_fun_calib_ntd.R"))
theme_set(theme_minimal())

opt_params_notimedepend <- readRDS(here("01_data/map_calib_ntd.RDS"))

Used overall deaths, prevalence of rx opioid use, and total OAT as calibration targets

Code
mod_opt_map_ntd <- mod_calib_ntd(as.numeric(opt_params_notimedepend$par))

8.6 Prevalence of Rx opioids use

Code
# prevalence of people on prescribed opioids

prop_opioids_rx_calib <- prop_opioids_rx_uncalib <- data.frame(matrix(NA, byrow = T,
                                             nrow = 12,
                                             ncol = length(2015:2018),
                                             list(NULL,
                                                  paste0("prop_opioids_rx_",
                                                         str_sub(2015:2018, 3, 4)))))
tot_pop_calib_tbl <- tot_pop_uncalib_tbl <- data.frame(matrix(NA, byrow = T,
                                         nrow = 12,
                                         ncol = length(2015:2018),
                                 list(NULL,
                                      paste0("tot_pop_",
                                             str_sub(2015:2018, 3, 4)))))

for (i in 1:length(2015:2018)) {
  yr <- 2014 + i
  prop_opioids_rx_uncalib[, i] <- assign(
    paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")]) /
      rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
  )
  
  tot_pop_uncalib_tbl[, i] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  
  prop_opioids_rx_calib[, i] <- assign(
    paste0("prop_opioids_rx_", str_sub(yr, 3, 4)),
    rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")]) /
      rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
  )
  
  tot_pop_calib_tbl[, i] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  }

prop_opioids_rx_uncalib_tbl <- prop_opioids_rx_uncalib %>% 
  pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>% 
  arrange(grp) %>% 
  mutate(year = rep(2015:2018,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_uncalib_tbl %>% 
              pivot_longer(1:4, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop))


prop_opioids_rx_calib_tbl <- prop_opioids_rx_calib %>% 
  pivot_longer(1:4, names_to = "grp", values_to = "target_val") %>% 
  arrange(grp) %>% 
  mutate(year = rep(2015:2018,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_calib_tbl %>% 
              pivot_longer(1:4, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop))


prop_opioids_rx_target_tbl <- calib_target_tbl %>% 
              filter(group == "prev_on_opioidsrx") %>% 
              select(year, target, group) %>% 
              rename(grp = group,
                     target_val = target) %>% 
              mutate(year_mon = paste(year, month.abb[6], sep = "_"))

prop_opioids_rx_uncalib_wei_mean <- prop_opioids_rx_uncalib_tbl %>% 
  group_by(year) %>% 
  summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
  ungroup() %>% 
  mutate(grp = "Uncalibrated Model") %>% 
  bind_rows(., prop_opioids_rx_target_tbl %>% 
              select(year, target_val) %>% 
              mutate(grp = "target")) %>% 
  bind_rows(., prop_opioids_rx_calib_tbl %>% 
              group_by(year) %>% 
              summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
              ungroup() %>% 
              mutate(grp = "Calibrated Model"))

prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))
prop_opioids_rx_uncalib_tbl$year_mon <- factor(prop_opioids_rx_uncalib_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))
prop_opioids_rx_calib_tbl$year_mon <- factor(prop_opioids_rx_calib_tbl$year_mon,
                                              levels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(2015:2018,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))

prop_opioids_rx_tbl <- prop_opioids_rx_uncalib_tbl %>% 
  mutate(group = "Uncalibrated Model") %>% 
  bind_rows(.,prop_opioids_rx_calib_tbl %>% 
              mutate(group = "Calibrated Model"))

wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>% 
  mutate(grp = "Target") %>% 
  bind_rows(., prop_opioids_rx_uncalib_tbl %>% 
               group_by(year) %>% 
               summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
                      grp = "Uncalibrated Model")) %>% 
  bind_rows(., prop_opioids_rx_calib_tbl %>% 
               group_by(year) %>% 
               summarise(target_val = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(year_mon = prop_opioids_rx_target_tbl$year_mon,
                      grp = "Calibrated Model"))

prop_opioids_rx_tbl$group <- factor(prop_opioids_rx_tbl$group,
                                  labels = c("Calibrated Model", "Uncalibrated Model"),
                                  levels = c("Calibrated Model", "Uncalibrated Model"))


p1 <- ggplot() +
  geom_point(data = prop_opioids_rx_tbl, aes(x = year_mon, y = target_val,
                                             color = group, fill = group)) +
  geom_point(data = wprop_opioids_rx_tbl %>% filter(grp == "Target"),
             aes(x = year_mon, y = target_val, color = grp, fill = grp), size = 2) +
   scale_color_brewer(palette = "Dark2") +
  scale_fill_brewer(palette = "Dark2") +
  xlab("Year_Month") + ylab("Prevalence of Rx opioid use") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

plotly::ggplotly(p1)

8.6.0.1 Comparison between uncalibrated annual weighted average and observed targets of prevalence of Rx opioids use

Code
wprop_opioids_rx_tbl$grp <- factor(wprop_opioids_rx_tbl$grp,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p2 <- ggplot() +
 geom_point(data = wprop_opioids_rx_tbl,
             aes(x = year_mon, y = target_val,
                 color = grp, fill = factor(grp)), size = 2) +
  xlab("Year") + ylab("Prevalence of Rx opioid use") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
  # ylim(0,0.2) +

plotly::ggplotly(p2)

8.7 Deaths and OD-Deaths

Code
################ Deaths 
# number of deaths

num_deaths_calib <- num_deaths_uncalib <- rep(NA, length(2016:2020))
for (i in 1:length(2016:2020)){
  yr <- 2015 + i
  num_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                       year_mon_cycle_tbl$mon == 12] + 1,
                                            "BO_DEATH"] -
    mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_DEATH"]
  
  num_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                       year_mon_cycle_tbl$mon == 12] + 1,
                                            "BO_DEATH"] -
    mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_DEATH"]
  }
  
# number of OD-deaths
  
num_od_deaths_calib <- num_od_deaths_uncalib <- rep(NA, length(2016:2021))
for (i in 1:length(2016:2021)){
  yr <- 2015 + i 
  num_od_deaths_uncalib[i] <- mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  num_od_deaths_calib[i] <- mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
}

deaths_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% c("total_deaths", "total_od_deaths")) %>% 
  mutate(target = ifelse(target == "total_deaths", "Total deaths",
                         ifelse(target == "total_od_deaths",
                                "Total OD-deaths", NA)),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_deaths_uncalib) %>% 
              mutate(target = "Total deaths") %>% 
              bind_rows(., data.frame(target_val = num_od_deaths_uncalib) %>%
                          mutate(target = "Total OD-deaths")) %>% 
              mutate(year = c(2016:2020, 2016:2021),
                     group = "Uncalibrated Model")) %>% 
   bind_rows(., data.frame(target_val = num_deaths_calib) %>% 
              mutate(target = "Total deaths") %>% 
              bind_rows(., data.frame(target_val = num_od_deaths_calib) %>%
                          mutate(target = "Total OD-deaths")) %>% 
              mutate(year = c(2016:2020, 2016:2021),
                     group = "Calibrated Model"))
  
deaths_target_uncalib_tbl$group <- factor(deaths_target_uncalib_tbl$group,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p3 <- ggplot() +
  geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Deaths") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
# +
#   facet_wrap(~target, scales = "free") + theme(legend.position="bottom")



plotly::ggplotly(p3)
Code
p4 <- ggplot() +
  geom_point(data = deaths_target_uncalib_tbl %>% filter(target == "Total OD-deaths"),
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total Opioid-related Overdose Deaths") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")


plotly::ggplotly(p4)

8.8 OAT

Code
# number of oat
  
num_oat_calib <- num_oat_uncalib <- rep(NA, length(2018:2021))
for (i in 1:length(2018:2021)){
  yr <- 2017 + i 
    num_oat_uncalib[i] <- sum(mod_basecase_notimedep$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                        year_mon_cycle_tbl$mon == 6] + 1,
                             c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
    
    num_oat_calib[i] <- sum(mod_opt_map_ntd$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                        year_mon_cycle_tbl$mon == 6] + 1,
                             c("BS_OAT_INI", "BS_OAT_MAINT", "BR_OAT_INI", "BR_OAT_MAINT")])
}


oat_target_uncalib_tbl <- calib_target_tbl %>% 
  select(year, target, group) %>% 
  rename(target_val = target,
         target = group) %>% 
  filter(target %in% c("total_oat")) %>% 
  mutate(target = ifelse(target == "total_oat",
                                "Total OAT", NA),
         group = "Target") %>% 
  bind_rows(., data.frame(target_val = num_oat_uncalib) %>%
                          mutate(target = "Total OAT",
                                 year = c(2018:2021),
                                 group = "Uncalibrated Model")) %>% 
   bind_rows(., data.frame(target_val = num_oat_calib) %>%
                          mutate(target = "Total OAT",
                                 year = c(2018:2021),
                                 group = "Calibrated Model"))

oat_target_uncalib_tbl$group <- factor(oat_target_uncalib_tbl$group,
                                   labels = c("Target", "Calibrated Model", "Uncalibrated Model"),
                                   levels = c("Target", "Calibrated Model", "Uncalibrated Model"))

p5 <- ggplot() +
  geom_point(data = oat_target_uncalib_tbl,
             aes(x = year, y = target_val, color = factor(group),
                 fill = factor(group))) +
  xlab("Year") + ylab("Total OAT counts") +
  scale_color_brewer(palette = "Dark2", name = "") +
  scale_fill_brewer(palette = "Dark2", name = "")
# +
#   facet_wrap(~target, scales = "free") + theme(legend.position="bottom")



plotly::ggplotly(p5)
Code
cbind.data.frame(parameter = calib_param_tbl_ntd$var_params,
                 basevalue = calib_param_tbl_ntd$basevalue,
                 opt_map = opt_params_notimedepend$par) %>% 
  filter(basevalue != opt_map)
                      parameter  basevalue    opt_map
1          p__BN_PN__BPO_MISUSE 0.00032511 0.00220011
2      p__BPO_CHRONIC__BO_OD_RX 0.00097898 0.00160398
3   p__BPO_PALLIATIVE__BO_DEATH 0.37500000 0.38250000
4 p__BS_OAT_MAINT__BS_OAT_MAINT 0.70000000 0.77000000

9 Appendix 4: Comparison between SIR calibrated model, uncalibrated model and calibration targets

Code
library(here)
library(plotly)
source(here("02_scripts/01_fun_data.R"))
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("03_outputs/03_map_calibmod_vs_targets_td.R"))

calib_samples <- readRDS(here("01_data/calib_samples_sir.RDS"))
uncalib_samples <- readRDS(here("01_data/uncalib_sample_sir.RDS"))

opt_params_sir_prev <- readRDS(here("01_data/prev_outcome_data_sir.RDS"))
opt_params_sir_deaths <- readRDS(here("01_data/death_outcome_data_sir.RDS"))
opt_params_sir_oddeaths <- readRDS(here("01_data/oddeath_outcome_data_sir.RDS"))
opt_params_sir_oat <- readRDS(here("01_data/oat_outcome_data_sir.RDS"))

params_sir_prev_uncalib <- readRDS(here("01_data/prev_outcome_data_sir_uncalib.RDS"))
params_sir_deaths_uncalib <- readRDS(here("01_data/death_outcome_data_sir_uncalib.RDS"))
params_sir_oddeaths_uncalib <- readRDS(here("01_data/oddeath_outcome_data_sir_uncalib.RDS"))
params_sir_oat_uncalib <- readRDS(here("01_data/oat_sir_uncalib.RDS"))

SIR calibration (calib_mod): first i sampled from prior 10000 sample sets, then i used them to calculated outcomes and likelihoods for each sample set, then i resampled from those 10000 with replacement with likelihood weights 10000 samples (unique sets = 6217), and ran the model on all of them and generated distributions for prevalence of rx opioid use, deaths, od deaths, total oat —- likelihood included deaths (gamma distribution), prevalence of rx opioid use (normal distribution instead of binomial), and total OAT

9.1 Prevalence of Rx opioid use

Code
dt_prev <- as.data.frame(params_sir_prev_uncalib) %>% 
  # rename(yr15 = "2015",
  #        yr16 = "2016",
  #        yr17 = "2017",
  #        yr18 = "2018") %>% 
  pivot_longer(1:4, names_to = "prev_year", values_to = "value") %>% 
  mutate(group = "Uncalibrated Sample") %>% 
  bind_rows(., as.data.frame(opt_params_sir_prev) %>%
              # rename(yr15 = "2015",
              #        yr16 = "2016",
              #        yr17 = "2017",
              #        yr18 = "2018") %>% 
              pivot_longer(1:4, names_to = "prev_year", values_to = "value") %>% 
              mutate(group = "SIR Calibrated Sample"))

prev_year = c("yr15", "yr16", "yr17", "yr18")

wprop_opioids_rx_tbl <- prop_opioids_rx_target_tbl %>% 
  rename(value = target_val) %>% 
  mutate(grp = "Target",
         prev_year = ifelse(grepl("15", year_mon),2015,
                            ifelse(grepl("16", year_mon), 2016,
                                   ifelse(grepl("17", year_mon), 2017,
                                          ifelse(grepl("18", year_mon), 2018, year_mon)))),
         prev_year = factor(prev_year)) %>%
  select(-year_mon) %>% 
  bind_rows(., prop_opioids_rx_uncalib_tbl %>% 
               group_by(year) %>% 
               summarise(value = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(prev_year = factor(year),
                      grp = "Uncalibrated Model - Point estimate")) %>% 
  bind_rows(., prop_opioids_rx_calib_tbl %>% 
               group_by(year) %>% 
               summarise(value = weighted.mean(target_val, tot_pop_val)) %>% 
               ungroup() %>% 
               mutate(prev_year = factor(year),
                      grp = "MAP Calibrated Model")) %>% 
  bind_rows(., dt_prev %>% 
              group_by(prev_year) %>% 
              summarize(value = mean(value)) %>% 
              mutate(grp = "SIR Calibrated Model - Mean",
                     prev_year = factor(prev_year))) %>% 
  bind_rows(., dt_prev %>% 
              group_by(prev_year) %>% 
              summarize(value = median(value)) %>% 
              mutate(grp = "SIR Calibrated Model - Median",
                     prev_year = factor(prev_year)))


ggplotly(ggplot(data = dt_prev,
       aes(y = value, x = prev_year)) +
  geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
  geom_point(data = wprop_opioids_rx_tbl,
             aes(y = value, x = prev_year, color = grp))+
  xlab("Year") +
  scale_color_brewer(palette = "Set1", name = "") +
  scale_fill_brewer(palette = "Set1", name = "") +
  facet_wrap(~group))
Code
ggplotly(ggplot(data = dt_prev,
                        aes(x = value)) +
                   geom_histogram(aes(fill = group), color="#e9ecef",
                                  alpha=0.3, position = 'identity')  +
                   geom_vline(data = wprop_opioids_rx_tbl, aes(xintercept = value, color = grp))+
                   theme(axis.text.x = element_text(angle = 45)) +
                   xlab("Prevalence of rx opioids use") + ylab("") +
                   scale_color_brewer(palette = "Set1", name = "") +
                   scale_fill_brewer(palette = "Set1", name = "") +
                   facet_wrap(~prev_year, scales = "free"))
Code
p_prev_sir <- ggplot(data = dt_prev,
       aes(y = value, x = prev_year)) +
  geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
  geom_point(data = wprop_opioids_rx_tbl,
             aes(y = value, x = prev_year, color = grp))+
  xlab("Year") +
  scale_color_brewer(palette = "Set1", name = "") +
  scale_fill_brewer(palette = "Set1", name = "") +
  facet_wrap(~group)
ggsave(here("04_report/p_prev_sir.png"))

9.2 Total Deaths

Code
dt_deaths <- as.data.frame(params_sir_deaths_uncalib) %>% 
  pivot_longer(1:5, names_to = "death_year", values_to = "value") %>% 
  mutate(group = "Uncalibrated Sample") %>% 
  bind_rows(., as.data.frame(opt_params_sir_deaths) %>% 
              pivot_longer(1:5, names_to = "death_year", values_to = "value") %>% 
              mutate(group = "SIR Calibrated Sample"))

ovrall_deaths_target_uncalib_tbl <- deaths_target_uncalib_tbl %>% 
  filter(target == "Total deaths") %>% 
  rename(value = target_val) %>% 
  mutate(grp = ifelse(group == "Calibrated Model", "MAP Calibrated Model",
                      ifelse(group == "Uncalibrated Model", "Uncalibrated Model - Point estimate",
                             ifelse(group == "Target", "Target", group))),
         death_year = factor(year)) %>% 
  select(-c(group, year)) %>% 
  bind_rows(., dt_deaths %>% 
              group_by(death_year) %>% 
              summarize(value = mean(value)) %>% 
              mutate(grp = "SIR Calibrated Model - Mean",
                     death_year = factor(death_year))) %>% 
  bind_rows(., dt_deaths %>% 
              group_by(death_year) %>% 
              summarize(value = median(value)) %>% 
              mutate(grp = "SIR Calibrated Model - Median",
                     death_year = factor(death_year)))


ggplotly(ggplot(data = dt_deaths, 
                aes(y = value, x = factor(death_year))) + 
           geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
           geom_point(data = ovrall_deaths_target_uncalib_tbl,
                      aes(y = value, x = factor(death_year), color = grp))+
           xlab("Year") +
           scale_color_brewer(palette = "Set1", name = "") +
           scale_fill_brewer(palette = "Set1", name = "") +
           facet_wrap(~group))
Code
ggplotly(ggplot(data = dt_deaths,
                        aes(x = value/1000)) +
                   geom_histogram(aes(fill = group), color="#e9ecef",
                                  alpha=0.4, position = 'identity')  +
                   # scale_fill_manual(values=c("#69b3a2", "#404080")) +
                   geom_vline(data = ovrall_deaths_target_uncalib_tbl,
                              aes(xintercept = value/1000, color = grp))+
                   theme(axis.text.x = element_text(angle = 45)) +
                   xlab("Total deaths (in thousands)")+ ylab("") +
           scale_color_brewer(palette = "Set1", name = "") +
                   scale_fill_brewer(palette = "Set1", name = "") +
                   xlim(NA, 320)+
                   facet_wrap(~death_year, scales = "free"))
Code
p_deaths_sir <- ggplot(data = dt_deaths, 
                aes(y = value, x = factor(death_year))) + 
           geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
           geom_point(data = ovrall_deaths_target_uncalib_tbl,
                      aes(y = value, x = factor(death_year), color = grp))+
           xlab("Year") +
           scale_color_brewer(palette = "Set1", name = "") +
           scale_fill_brewer(palette = "Set1", name = "") +
           facet_wrap(~group)
ggsave(here("04_report/p_deaths_sir.png"))

9.4 Total OAT counts

Code
dt_oat <- as.data.frame(params_sir_oat_uncalib) %>% 
  pivot_longer(1:length(2018:2021), names_to = "oat_year", values_to = "value") %>% 
  mutate(group = "Uncalibrated Sample") %>% 
  bind_rows(., as.data.frame(opt_params_sir_oat) %>% 
              pivot_longer(1:length(2018:2021), names_to = "oat_year", values_to = "value") %>% 
              mutate(group = "SIR Calibrated Sample"))

oat_target_uncalib_tbl <- oat_target_uncalib_tbl %>% 
  rename(value = target_val) %>% 
  mutate(grp = ifelse(group == "Calibrated Model", "MAP Calibrated Model",
                      ifelse(group == "Uncalibrated Model", "Uncalibrated Model - Point estimate",
                             ifelse(group == "Target", "Target", group))),
         oat_year = factor(year)) %>% 
  select(-c(group, year)) %>% 
  bind_rows(., dt_oat %>% 
              group_by(oat_year) %>% 
              summarize(value = mean(value)) %>% 
              mutate(grp = "SIR Calibrated Model - Mean",
                     oat_year = factor(oat_year))) %>% 
  bind_rows(., dt_oat %>% 
              group_by(oat_year) %>% 
              summarize(value = median(value)) %>% 
              mutate(grp = "SIR Calibrated Model - Median",
                     oat_year = factor(oat_year)))



ggplotly(ggplot(data = dt_oat, 
                aes(y = value, x = factor(oat_year))) + 
           geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
           geom_point(data = oat_target_uncalib_tbl,
                      aes(y = value, x = factor(oat_year), color = grp))+
           xlab("Year") +
           scale_color_brewer(palette = "Set1", name = "") +
           scale_fill_brewer(palette = "Set1", name = "") +
           facet_wrap(~group))
Code
ggplotly(ggplot(data = dt_oat,
                        aes(x = value)) +
                   geom_histogram(aes(fill = group), color="#e9ecef",
                                  alpha=0.4, position = 'identity')  +
                   geom_vline(data = oat_target_uncalib_tbl,
                              aes(xintercept = value, color = grp))+
                   theme(axis.text.x = element_text(angle = 45)) +
                   xlab("Total OAT")+ ylab("") +
           scale_color_brewer(palette = "Set1", name = "") +
                   scale_fill_brewer(palette = "Set1", name = "") +
                   facet_wrap(~oat_year, scales = "free"))
Code
p_oat_sir <- ggplot(data = dt_oat, 
                aes(y = value, x = factor(oat_year))) + 
           geom_violin(aes(fill = group), position = "dodge", alpha = 0.5) +
           geom_point(data = oat_target_uncalib_tbl,
                      aes(y = value, x = factor(oat_year), color = grp))+
           xlab("Year") +
           scale_color_brewer(palette = "Set1", name = "") +
           scale_fill_brewer(palette = "Set1", name = "") +
           facet_wrap(~group)

ggsave(here("04_report/p_oat_sir.png"))

10 Appendix 5: CBA results - point estimates

Code
library(here)
library(flextable)
library(plotly)
source(here("02_scripts/03a_fun_calib_td.R"))
source(here("02_scripts/06_interventions_models.R"))

theme_set(theme_minimal())

Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 Scenario 2: Scenario 2: Doesn’t account for Fentanyl contamination and has the same probability from Illicit use to OD and from OD to OD death

Code
trace_tbl <- as_tibble(mod_basecase$m_M) %>% 
  mutate(cycle_num = 0:n_cycles) %>% 
  pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
  mutate(mod = "No Interventions") %>% 
  bind_rows(., as_tibble(mod_nalx_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Naloxone")) %>% 
  bind_rows(., as_tibble(mod_ss_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Safer Supply")) %>% 
  bind_rows(., as_tibble(mod_pres_guid_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Prescription Guidelines")) %>% 
  bind_rows(., as_tibble(mod_all_int_1$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "All Interventions"))

trace_tbl$mod <- factor(trace_tbl$mod,
                        levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
                        labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))

v_extra_sevbi <- c(mod_basecase$extra_sev_bi, mod_nalx_1$extra_sev_bi,
                   mod_ss_1$extra_sev_bi, mod_pres_guid_1$extra_sev_bi,
                   mod_all_int_1$extra_sev_bi)

v_extra_modbi <- c(mod_basecase$extra_mod_bi, mod_nalx_1$extra_mod_bi,
                   mod_ss_1$extra_mod_bi, mod_pres_guid_1$extra_mod_bi,
                   mod_all_int_1$extra_mod_bi)

v_extra_od <- c(mod_basecase$extra_od_deaths, mod_nalx_1$extra_od_deaths,
                   mod_ss_1$extra_od_deaths, mod_pres_guid_1$extra_od_deaths,
                   mod_all_int_1$extra_od_deaths)

mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)

for (j in 1:length(mod_names)){
  
  trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
                      trace_tbl$cycle_num == 180 &
                      trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_OD_DEATH" &
                      trace_tbl$cycle_num == 180 &
                      trace_tbl$mod == mod_names[j]] + v_extra_od[j]
  
  trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_MOD_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] + v_extra_modbi[j]
  
  trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] <- trace_tbl$count[trace_tbl$state == "BO_SEVERE_BI" &
                      trace_tbl$cycle_num == 180  &
                      trace_tbl$mod == mod_names[j]] + v_extra_sevbi[j]
}

trace_tbl$state <- factor(trace_tbl$state,
                              levels = v_state_names,
                              labels = v_state_names)
Code
trace_tbl_2 <- as_tibble(mod_basecase_2$m_M) %>% 
  mutate(cycle_num = 0:n_cycles) %>% 
  pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
  mutate(mod = "No Interventions") %>% 
  bind_rows(., as_tibble(mod_nalx_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Naloxone")) %>% 
  bind_rows(., as_tibble(mod_ss_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Safer Supply")) %>% 
  bind_rows(., as_tibble(mod_pres_guid_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "Prescription Guidelines")) %>% 
  bind_rows(., as_tibble(mod_all_int_2$m_M) %>% 
              mutate(cycle_num = 0:n_cycles) %>% 
              pivot_longer(1:n_states, names_to = "state", values_to = "count") %>% 
              mutate(mod = "All Interventions"))

trace_tbl_2$mod <- factor(trace_tbl_2$mod,
                        levels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"),
                        labels = c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions"))

v_extra_sevbi_2 <- c(mod_basecase_2$extra_sev_bi, mod_nalx_2$extra_sev_bi,
                   mod_ss_2$extra_sev_bi, mod_pres_guid_2$extra_sev_bi,
                   mod_all_int_2$extra_sev_bi)

v_extra_modbi_2 <- c(mod_basecase_2$extra_mod_bi, mod_nalx_2$extra_mod_bi,
                   mod_ss_2$extra_mod_bi, mod_pres_guid_2$extra_mod_bi,
                   mod_all_int_2$extra_mod_bi)

v_extra_od_2 <- c(mod_basecase_2$extra_od_deaths, mod_nalx_2$extra_od_deaths,
                   mod_ss_2$extra_od_deaths, mod_pres_guid_2$extra_od_deaths,
                   mod_all_int_2$extra_od_deaths)

for (j in 1:length(mod_names)){
  
  trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
                      trace_tbl_2$cycle_num == 180 &
                      trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_OD_DEATH" &
                      trace_tbl_2$cycle_num == 180 &
                      trace_tbl_2$mod == mod_names[j]] + v_extra_od_2[j]
  
  trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_MOD_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] + v_extra_modbi_2[j]

    trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] <- trace_tbl_2$count[trace_tbl_2$state == "BO_SEVERE_BI" &
                      trace_tbl_2$cycle_num == 180  &
                      trace_tbl_2$mod == mod_names[j]] + v_extra_sevbi_2[j]
}

trace_tbl_2$state <- factor(trace_tbl_2$state,
                              levels = v_state_names,
                              labels = v_state_names)

trace_tbl <- trace_tbl %>% 
  mutate(scenario = "Scenario 1") %>% 
  bind_rows(., trace_tbl_2 %>% 
              mutate(scenario = "Scenario 2"))

10.1 Cumulative OD Deaths

10.1.0.0.1 Cumulative OD Deaths - Scenario 1
Code
p1 <- ggplot(trace_tbl %>% 
               filter(state == "BO_OD_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 1") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
               geom_line() +
               geom_point(size = 1) +
               ylab("Cumulative OD Deaths") + xlab("Year") +
               labs(colour = "Intervention")

ggplotly(p1)
10.1.0.0.2 Cumulative OD Deaths - Scenario 2
Code
p2 <- ggplot(trace_tbl %>% 
               filter(state == "BO_OD_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 2") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
               geom_line() +
               geom_point(size = 1) +
               ylab("Cumulative OD Deaths") + xlab("Year") +
               labs(colour = "Intervention")

ggplotly(p2)
10.1.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
               filter(state == "BO_OD_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               mutate(year = rep(c(2015:2029), length(mod_names)*2)),
             aes(x = year, y = count, colour = mod)) +
               geom_line() +
               geom_point(size = 1) +
               ylab("Cumulative OD Deaths") + xlab("Year") +
               labs(colour = "Intervention") +
               facet_wrap(~scenario) +
               labs(caption = "Scenario 1: Accounts for the effect of Fentanyl contamination starting from 2018 \n
                    Scenario 2: Doesn't account for Fentanyl contamination and has the same probability from \n Illicit use to OD and from OD to OD death"))

10.2 New OD Deaths per year

Code
inci_od_deaths_basecase <- inci_od_deaths_nalox <- inci_od_deaths_ss <- inci_od_deaths_pg <- inci_od_deaths_all <- inci_od_deaths_basecase_2 <- inci_od_deaths_nalox_2 <- inci_od_deaths_ss_2 <- inci_od_deaths_pg_2 <- inci_od_deaths_all_2 <-  rep(NA, 14)

for (i in 1:14){
  yr <- 2015 + i 
  inci_od_deaths_basecase[i] <- mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_basecase_2[i] <- mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_basecase_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_nalox[i] <- mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_nalx_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
   inci_od_deaths_nalox_2[i] <- mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_nalx_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_ss[i] <- mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_ss_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
   inci_od_deaths_ss_2[i] <- mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_ss_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_pg[i] <- mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_pres_guid_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
   inci_od_deaths_pg_2[i] <- mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_pres_guid_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  
  inci_od_deaths_all[i] <- mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_all_int_1$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
  
  inci_od_deaths_all_2[i] <- mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr &
                                                                          year_mon_cycle_tbl$mon == 12] + 1,
                                               "BO_OD_DEATH"] -
    mod_all_int_2$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr - 1 &
                                                year_mon_cycle_tbl$mon == 12] + 1,
                     "BO_OD_DEATH"]
}

inc_od_death_tbl <- data.frame(intervention = c(rep(mod_names, each = 14)),
           year = c(rep(2016:2029, 5)),
           inci_od_deaths = c(inci_od_deaths_basecase,
                              inci_od_deaths_nalox,
                              inci_od_deaths_ss,
                              inci_od_deaths_pg,
                              inci_od_deaths_all),
           scenario = rep("Scenario 1", 14*5)) %>% 
  bind_rows(., data.frame(intervention = c(rep(mod_names, each = 14)),
           year = c(rep(2016:2029, 5)),
           inci_od_deaths = c(inci_od_deaths_basecase_2,
                              inci_od_deaths_nalox_2,
                              inci_od_deaths_ss_2, 
                              inci_od_deaths_pg_2,
                              inci_od_deaths_all_2),
           scenario = rep("Scenario 2", 14*5)))

saveRDS(inc_od_death_tbl, file = here("01_data/inc_od_death_tbl_mod.RDS"))

target_od_deaths <- calib_target_tbl %>% 
  filter(group == "total_od_deaths") %>% 
  select(year, target) %>% 
  rename(inci_od_deaths = target) %>% 
  mutate(intervention = "Target")

saveRDS(target_od_deaths, file = here("01_data/inc_od_death_tbl_target.RDS"))

inc_od_death_tbl$intervention <- factor(inc_od_death_tbl$intervention,
                                        levels = mod_names, labels = mod_names)
10.2.0.0.1 Scenario 1
Code
p3 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 1"),
                        aes(x = year, y = inci_od_deaths, color = intervention)) +
               geom_line() +
               geom_point(size = 1) +
                 geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
               ylab("New OD Deaths per year") + xlab("Year") +
               labs(colour = "Intervention") +
                 scale_color_brewer(palette = "Dark2")

ggplotly(p3)
10.2.0.0.2 Scenario 2
Code
p4 <- ggplot(inc_od_death_tbl %>% filter(scenario == "Scenario 2"),
                        aes(x = year, y = inci_od_deaths, color = intervention)) +
               geom_line() +
               geom_point(size = 1) +
                 geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
               ylab("New OD Deaths per year") + xlab("Year") +
               labs(colour = "Intervention") +
                 scale_color_brewer(palette = "Dark2")

ggplotly(p4)
10.2.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(inc_od_death_tbl,
                        aes(x = year, y = inci_od_deaths, color = intervention)) +
               geom_line() +
               geom_point(size = 1) +
                 geom_point(data = target_od_deaths, aes (x = year, y = inci_od_deaths)) +
               ylab("New OD Deaths per year") + xlab("Year") +
               labs(colour = "Intervention") +
                 scale_color_brewer(palette = "Dark2") +
                 facet_wrap(~scenario))

10.3 Cumulative Deaths

10.3.0.0.1 Scenario 1
Code
p5 <- ggplot(trace_tbl %>% 
               filter(state == "BO_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 1") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
  geom_line() +
  ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
    labs(colour = "Intervention") 

ggplotly(p5)
10.3.0.0.2 Scenario 2
Code
p6 <- ggplot(trace_tbl %>% 
               filter(state == "BO_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               filter(scenario == "Scenario 2") %>% 
               mutate(year = rep(c(2015:2029), length(mod_names))),
             aes(x = year, y = count, colour = mod)) +
  geom_line() +
  ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
    labs(colour = "Intervention") 

ggplotly(p6)
10.3.0.0.3 Both scenarios
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
               filter(state == "BO_DEATH") %>% 
               filter(cycle_num != 0) %>% 
               filter((cycle_num + 1) %% 12 == 1) %>% 
               mutate(year = rep(c(2015:2029), length(mod_names)*2)),
             aes(x = year, y = count, colour = mod)) +
  geom_line() +
  ylab("Cumulative Deaths (excluding OD-deaths)") + xlab("Year") +
    labs(colour = "Intervention") +
    facet_wrap(~scenario))

10.4 Trace Plots

Code
set.seed(122345)
colors_pal <- Polychrome::createPalette(n_states,
                                        c("#084C61", "#DB504A",
                                                   "#E3B505", "#4F6D7A", 
                                                      "#56A3A6"))
names(colors_pal) <- NULL

ggplotly(ggplot(trace_tbl %>% 
         filter(grepl("BPO", state)), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% filter(cycle_num == 180) %>% 
         filter(grepl("BPO", state)), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal) +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BPO_MISUSE", "BI_ILLICIT", 
                             "BS_OAT_INI", "BS_OAT_MAINT")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BPO_MISUSE", "BI_ILLICIT",
                             "BS_OAT_INI", "BS_OAT_MAINT")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal)  +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BS_OAT_SS", "BS_SS")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BS_OAT_SS", "BS_SS")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal) +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
                             "BO_MOD_BI", "BO_SEVERE_BI")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BO_OD_ILLICIT", "BO_OD_RX",
                             "BO_MOD_BI", "BO_SEVERE_BI")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal)  +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))
Code
plotly::ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
                             "BR_OAT_MAINT", "BR_OAT_SS",
                             "BR_SS", "BR_OD_ILLICIT")), aes(x = cycle_num, y = count, colour = state)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  ggrepel::geom_text_repel(data = trace_tbl %>% 
         filter(state %in% c("BR_ILLICIT", "BR_OAT_INI",
                             "BR_OAT_MAINT", "BR_OAT_SS",
                             "BR_SS", "BR_OD_ILLICIT")) %>% filter(cycle_num == 180), 
                  aes(x = cycle_num, y = count, label = state), 
                  check_overlap = T, size = 3) +
  scale_color_manual(values = colors_pal)  +
  # facet_wrap(~mod) +
    facet_grid(rows = vars(mod), cols = vars(scenario),
               switch = "y"))

10.5 CBA results

10.5.0.1 Scenario 1

Code
cost_diff <- c(mod_nalx_1$total_net_present_cost - mod_basecase$total_net_present_cost,
               mod_ss_1$total_net_present_cost - mod_basecase$total_net_present_cost,
               mod_pres_guid_1$total_net_present_cost - mod_basecase$total_net_present_cost,
               mod_all_int_1$total_net_present_cost - mod_basecase$total_net_present_cost)

cost_diff_per <- c(round((cost_diff/c(mod_basecase$total_net_present_cost))*100,2))

deaths_eff_tbl <- trace_tbl %>% 
  filter(scenario == "Scenario 1") %>% 
  filter(cycle_num == 180) %>% 
  select(-cycle_num) %>% 
  pivot_wider(names_from = mod, values_from = count) %>% 
  filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>% 
  mutate(nalx_effect = Naloxone - `No Interventions`,
         ss_effect = `Safer Supply` - `No Interventions`,
         pg_effect = `Prescription Guidelines` - `No Interventions`,
         all_effect = `All Interventions` - `No Interventions`,
         nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
         ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
         pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
         all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>% 
  select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))


total_death_oddeaths <- trace_tbl %>% 
  filter(scenario == "Scenario 1") %>%
  filter(cycle_num == 180) %>% 
  select(-cycle_num) %>% 
  pivot_wider(names_from = mod, values_from = count) %>% 
  filter(state %in% c("BO_DEATH", "BO_OD_DEATH"))

od_deaths_diff <- c(deaths_eff_tbl$nalx_effect[1], deaths_eff_tbl$ss_effect[1],
                 deaths_eff_tbl$pg_effect[1], deaths_eff_tbl$all_effect[1])
od_deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[1], deaths_eff_tbl$ss_effect_per[1],
                 deaths_eff_tbl$pg_effect_per[1], deaths_eff_tbl$all_effect_per[1])
deaths_diff <- c(deaths_eff_tbl$nalx_effect[2], deaths_eff_tbl$ss_effect[2],
                 deaths_eff_tbl$pg_effect[2], deaths_eff_tbl$all_effect[2])
deaths_diff_per <- c(deaths_eff_tbl$nalx_effect_per[2], deaths_eff_tbl$ss_effect_per[2],
                 deaths_eff_tbl$pg_effect_per[2], deaths_eff_tbl$all_effect_per[2])


saveRDS(cbind.data.frame(cost_diff, cost_diff_per,
                         od_deaths_diff, od_deaths_diff_per,
                         deaths_diff, deaths_diff_per), here("01_data/point_estiamte.RDS"))
saveRDS(total_death_oddeaths, here("01_data/total_deaths_oddeaths.RDS"))

costs <- paste0(round(cost_diff/1000000, 0), " (", cost_diff_per, "%)")
deaths <- paste0(round(deaths_diff, 0), " (", deaths_diff_per, "%)")
oddeaths <- paste0(round(od_deaths_diff, 0), " (", od_deaths_diff_per, "%)")

effects_tbl <- cbind.data.frame(mod_names[-1], costs, deaths, oddeaths)

effects_tbl |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Mean Change Compared to No Intervention")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         costs = "Discounted Net Present Costs \n in Millions (%)",
                         deaths = "Deaths (%)",
                         oddeaths = "Overdose-related Deaths (%)")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Mean Change Compared to No Intervention

Interventions

Discounted Net Present Costs
in Millions (%)

Deaths (%)

Overdose-related Deaths (%)

Naloxone

102 (0.01%)

396 (0.01%)

-6420 (-4.05%)

Safer Supply

4233 (0.53%)

-13956 (-0.31%)

-3138 (-1.98%)

Prescription Guidelines

-26103 (-3.28%)

15252 (0.34%)

-1764 (-1.11%)

All Interventions

-24876 (-3.13%)

14115 (0.31%)

-8543 (-5.39%)

10.5.0.2 Scenario 2

Code
cost_diff_2 <- c(mod_nalx_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
               mod_ss_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
               mod_pres_guid_2$total_net_present_cost - mod_basecase_2$total_net_present_cost,
               mod_all_int_2$total_net_present_cost - mod_basecase_2$total_net_present_cost)

cost_diff_per_2 <- c(round((cost_diff_2/c(mod_basecase_2$total_net_present_cost))*100,2))

deaths_eff_tbl_2 <- trace_tbl %>% 
  filter(scenario == "Scenario 2") %>% 
  filter(cycle_num == 180) %>% 
  select(-cycle_num) %>% 
  pivot_wider(names_from = mod, values_from = count) %>% 
  filter(state %in% c("BO_DEATH", "BO_OD_DEATH")) %>% 
  mutate(nalx_effect = Naloxone - `No Interventions`,
         ss_effect = `Safer Supply` - `No Interventions`,
         pg_effect = `Prescription Guidelines` - `No Interventions`,
         all_effect = `All Interventions` - `No Interventions`,
         nalx_effect_per = round(((Naloxone - `No Interventions`)/`No Interventions`)*100, 2),
         ss_effect_per = round(((`Safer Supply` - `No Interventions`)/`No Interventions`)*100, 2),
         pg_effect_per = round(((`Prescription Guidelines` - `No Interventions`)/`No Interventions`)*100, 2),
         all_effect_per = round(((`All Interventions` - `No Interventions`)/`No Interventions`)*100, 2)) %>% 
  select(-c(Naloxone, `No Interventions`, `Safer Supply`, `Prescription Guidelines`, `All Interventions`))

od_deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[1], deaths_eff_tbl_2$ss_effect[1],
                 deaths_eff_tbl_2$pg_effect[1], deaths_eff_tbl_2$all_effect[1])
od_deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[1], deaths_eff_tbl_2$ss_effect_per[1],
                 deaths_eff_tbl_2$pg_effect_per[1], deaths_eff_tbl_2$all_effect_per[1])
deaths_diff_2 <- c(deaths_eff_tbl_2$nalx_effect[2], deaths_eff_tbl_2$ss_effect[2],
                 deaths_eff_tbl_2$pg_effect[2], deaths_eff_tbl_2$all_effect[2])
deaths_diff_per_2 <- c(deaths_eff_tbl_2$nalx_effect_per[2], deaths_eff_tbl_2$ss_effect_per[2],
                 deaths_eff_tbl_2$pg_effect_per[2], deaths_eff_tbl_2$all_effect_per[2])

costs_2 <- paste0(round(cost_diff_2/1000000, 0), " (", cost_diff_per_2, "%)")
deaths_2 <- paste0(round(deaths_diff_2, 0), " (", deaths_diff_per_2, "%)")
oddeaths_2 <- paste0(round(od_deaths_diff_2, 0), " (", od_deaths_diff_per_2, "%)")

effects_tbl_2 <- cbind.data.frame(mod_names[-1], costs_2, deaths_2, oddeaths_2)

effects_tbl_2 |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Mean Change Compared to No Intervention")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         costs_2 = "Discounted Net Present Costs \n in Millions (%)",
                         deaths_2 = "Deaths (%)",
                         oddeaths_2 = "Overdose-related Deaths (%)")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Mean Change Compared to No Intervention

Interventions

Discounted Net Present Costs
in Millions (%)

Deaths (%)

Overdose-related Deaths (%)

Naloxone

85 (0.01%)

303 (0.01%)

-4861 (-3.9%)

Safer Supply

4047 (0.51%)

-12743 (-0.28%)

-2435 (-1.95%)

Prescription Guidelines

-26098 (-3.29%)

15308 (0.34%)

-1485 (-1.19%)

All Interventions

-24882 (-3.13%)

14105 (0.31%)

-6607 (-5.3%)

Code
tbl_effect_plot1 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
                                    deaths_diff_per, od_deaths_diff_per) %>% 
  pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>% 
  mutate(scenario = "Scenario 1") %>% 
  bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
                                    deaths_diff_per_2, od_deaths_diff_per_2) %>% 
  pivot_longer(2:4, names_to = "group", values_to = "percent_change") %>% 
  mutate(scenario = "Scenario 2")) %>% 
  mutate(group = ifelse(group == "cost_diff_per_2", "cost_diff_per",
                        ifelse(group == "deaths_diff_per_2", "deaths_diff_per",
                               ifelse(group == "od_deaths_diff_per_2", "od_deaths_diff_per", group))))
10.5.0.2.1 Scenario 1
Code
p7 <- ggplot(data = tbl_effect_plot1 %>% 
               filter(scenario == "Scenario 1"), aes(x = group, y = percent_change))+
  geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","OD-releated Deaths"))+
  theme_minimal()

ggplotly(p7)
10.5.0.2.2 Scenario 2
Code
p8 <- ggplot(data = tbl_effect_plot1 %>% 
               filter(scenario == "Scenario 2"), aes(x = group, y = percent_change))+
  geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","OD-releated Deaths"))+
  theme_minimal()

ggplotly(p8)
10.5.0.2.3 Both
Code
ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = percent_change))+
  geom_point(aes(colour = Intervention, shape = scenario), size = 1.5) +
  scale_color_brewer(palette = "Dark2") +
  # scale_shape_manual(values = c(15, 16, 17, 18)) +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","OD-releated Deaths"))+
    theme_bw() +
    facet_wrap(~scenario))
Code
tbl_effect_plot2 <- cbind.data.frame(Intervention = mod_names[-1], cost_diff_per,
                                     od_deaths_diff_per) %>% 
  mutate(scenario = "Scenario 1") %>% 
  bind_rows(., cbind.data.frame(Intervention = mod_names[-1], cost_diff_per_2,
                                     od_deaths_diff_per_2) %>%
              rename(cost_diff_per = "cost_diff_per_2",
                     od_deaths_diff_per = "od_deaths_diff_per_2") %>% 
              mutate(scenario = "Scenario 2"))
10.5.0.2.4 Scenario 1
Code
p9 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 1"),
             aes(x = cost_diff_per, y = od_deaths_diff_per))+
  geom_point(aes(color = Intervention, shape = scenario)) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  xlab("Mean Change in Costs Compared to No Intervention (%)") +
  ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
  theme_minimal()

ggplotly(p9)
10.5.0.2.5 Scenario 2
Code
p10 <- ggplot(data = tbl_effect_plot2 %>% filter(scenario == "Scenario 2"),
             aes(x = cost_diff_per, y = od_deaths_diff_per))+
  geom_point(aes(color = Intervention, shape = scenario)) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  xlab("Mean Change in Costs Compared to No Intervention (%)") +
  ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
  theme_minimal()

ggplotly(p10)
10.5.0.2.6 Both
Code
ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_diff_per, y = od_deaths_diff_per))+
  geom_point(aes(color = Intervention, shape = scenario)) +
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  geom_vline(xintercept = 0) +
  xlab("Mean Change in Costs Compared to No Intervention (%)") +
  ylab("Mean Change in Overdose-related deaths \n Compared to No Intervention (%)") +
  theme_bw() +
    facet_wrap(~scenario))

10.6 Secondary Outcomes

Code
v_yrs_prev <- c(2015:2029)
yrs_prev <- length(v_yrs_prev)


prev_fun <- function(mod, i, v_states){
  
  prop_uncalib <- data.frame(matrix(NA,byrow = T, nrow = 12, 
                                             ncol = yrs_prev,
                                             list(NULL,
                                                  paste0("prop_",
                                                         str_sub(v_yrs_prev, 3, 4)))))
  
  tot_pop_uncalib_tbl <- data.frame(matrix(NA,byrow = T, nrow = 12, 
                                             ncol = yrs_prev,
                                 list(NULL,
                                      paste0("tot_pop_",
                                             str_sub(v_yrs_prev, 3, 4)))))
  
  for (j in 1:yrs_prev) {
    yr <- (v_yrs_prev[1] - 1) + j
    
    
    if (length(v_states) > 1)
      {
      prop_uncalib[, j] <- assign(
        paste0("prop_", str_sub(yr, 3, 4)),
        rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             v_states]) /
          rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
        )
    } else {
      prop_uncalib[, j] <- assign(
      paste0("prop_", str_sub(yr, 3, 4)),
      mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,
                             v_states]) /
      rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
    
    }
    
    
  
  tot_pop_uncalib_tbl[, j] <- assign(
    paste0("tot_pop_", str_sub(yr, 3, 4)),
    rowSums(mod$m_M[year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year == yr] + 1,])
  )
  
  }
  
  # final table for monthly prevalence over 15 years
  prop_uncalib_tbl <- prop_uncalib %>% 
    pivot_longer(1:15, names_to = "grp", values_to = "value") %>% 
    arrange(grp) %>% 
    mutate(year = rep(v_yrs_prev,
                              each = 12),
         year_mon = paste(year, month.abb[1:12], sep = "_")) %>% 
  bind_cols(., tot_pop_uncalib_tbl %>% 
              pivot_longer(1:15, names_to = "tot_pop",
                           values_to = "tot_pop_val") %>% 
              arrange(tot_pop) %>% select(-tot_pop)) %>% 
    mutate(intervention = mod_names[i])
  
  # table with weighted yearly prevalence
  prop_uncalib_wei_mean <- prop_uncalib_tbl %>% 
    group_by(year) %>% 
    summarise(value = weighted.mean(value, tot_pop_val)) %>% 
    ungroup() %>% 
    mutate(intervention = mod_names[i]) 
  
  return (list(prop_uncalib_tbl = prop_uncalib_tbl,
               prop_uncalib_wei_mean = prop_uncalib_wei_mean))
}

10.6.1 Severe Brain injury

Code
noint_prev_sevbi_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
noint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
nalx_prev_sevbi_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
nalx_prev_sevbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
ss_prev_sevbi_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
ss_prev_sevbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
pg_prev_sevbi_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
pg_prev_sevbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean
allint_prev_sevbi_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_tbl
allint_prev_sevbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = "BO_SEVERE_BI")$prop_uncalib_wei_mean


prev_sevbi_mon_tbl <- noint_prev_sevbi_mon_tbl %>% 
  bind_rows(., nalx_prev_sevbi_mon_tbl, ss_prev_sevbi_mon_tbl,
            pg_prev_sevbi_mon_tbl, allint_prev_sevbi_mon_tbl)

prev_sevbi_yr_tbl <- noint_prev_sevbi_yr_tbl %>% 
  bind_rows(., nalx_prev_sevbi_yr_tbl, ss_prev_sevbi_yr_tbl,
            pg_prev_sevbi_yr_tbl, allint_prev_sevbi_yr_tbl) 

prev_sevbi_mon_tbl$year_mon <- factor(prev_sevbi_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% "BO_SEVERE_BI") %>%
           filter(scenario == "Scenario 1") %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_sevbi_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of severe brain injury") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.2 Prevalence of moderate Brain injury

Code
noint_prev_modbi_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_MOD_BI")$prop_uncalib_tbl
noint_prev_modbi_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
nalx_prev_modbi_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = "BO_MOD_BI")$prop_uncalib_tbl
nalx_prev_modbi_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
ss_prev_modbi_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = "BO_MOD_BI")$prop_uncalib_tbl
ss_prev_modbi_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
pg_prev_modbi_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = "BO_MOD_BI")$prop_uncalib_tbl
pg_prev_modbi_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean
allint_prev_modbi_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = "BO_MOD_BI")$prop_uncalib_tbl
allint_prev_modbi_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = "BO_MOD_BI")$prop_uncalib_wei_mean


prev_modbi_mon_tbl <- noint_prev_modbi_mon_tbl %>% 
  bind_rows(., nalx_prev_modbi_mon_tbl, ss_prev_modbi_mon_tbl,
            pg_prev_modbi_mon_tbl, allint_prev_modbi_mon_tbl)

prev_modbi_yr_tbl <- noint_prev_modbi_yr_tbl %>% 
  bind_rows(., nalx_prev_modbi_yr_tbl, ss_prev_modbi_yr_tbl,
            pg_prev_modbi_yr_tbl, allint_prev_modbi_yr_tbl) 

prev_modbi_mon_tbl$year_mon <- factor(prev_modbi_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% "BO_MOD_BI") %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_modbi_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of moderate brain injury") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.3 Prevalence of substance use treatment

Code
v_states_prev_tx <- c("BS_DETOX", "BS_OAT_INI", "BS_OAT_MAINT",
                          "BS_OAT_SS", "BS_SS", "BR_OAT_INI", "BR_OAT_MAINT",
                          "BR_OAT_SS", "BR_SS")

noint_prev_tx_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_prev_tx)$prop_uncalib_tbl
noint_prev_tx_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
nalx_prev_tx_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = v_states_prev_tx)$prop_uncalib_tbl
nalx_prev_tx_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
ss_prev_tx_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = v_states_prev_tx)$prop_uncalib_tbl
ss_prev_tx_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
pg_prev_tx_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = v_states_prev_tx)$prop_uncalib_tbl
pg_prev_tx_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean
allint_prev_tx_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = v_states_prev_tx)$prop_uncalib_tbl
allint_prev_tx_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_prev_tx)$prop_uncalib_wei_mean


prev_tx_mon_tbl <- noint_prev_tx_mon_tbl %>% 
  bind_rows(., nalx_prev_tx_mon_tbl, ss_prev_tx_mon_tbl,
            pg_prev_tx_mon_tbl, allint_prev_tx_mon_tbl)

prev_tx_yr_tbl <- noint_prev_tx_yr_tbl %>% 
  bind_rows(., nalx_prev_tx_yr_tbl, ss_prev_tx_yr_tbl,
            pg_prev_tx_yr_tbl, allint_prev_tx_yr_tbl) 

prev_tx_mon_tbl$year_mon <- factor(prev_tx_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))


ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% v_states_prev_tx) %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_tx_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of substance use treatment") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.4 Prevalence of prescription opioid use

Code
prop_opioids_rx_target_tbl <- calib_target_tbl %>% 
              filter(group == "prev_on_opioidsrx") %>% 
              select(year, target) %>% 
              rename(value = target) %>% 
              mutate(year_mon = paste(year, month.abb[7], sep = "_"))
  
v_states_prev_rx_opd <- c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")


noint_prev_rx_opd_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
noint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
nalx_prev_rx_opd_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
nalx_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
ss_prev_rx_opd_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
ss_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
pg_prev_rx_opd_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
pg_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
allint_prev_rx_opd_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_tbl
allint_prev_rx_opd_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_prev_rx_opd)$prop_uncalib_wei_mean
  

prev_rx_opd_mon_tbl <- noint_prev_rx_opd_mon_tbl %>% 
  bind_rows(., nalx_prev_rx_opd_mon_tbl, ss_prev_rx_opd_mon_tbl,
            pg_prev_rx_opd_mon_tbl, allint_prev_rx_opd_mon_tbl) %>% 
    bind_rows(., prop_opioids_rx_target_tbl %>% 
              select(year, value, year_mon) %>% 
              mutate(intervention = "Target"))

prev_rx_opd_yr_tbl <- noint_prev_rx_opd_yr_tbl %>% 
  bind_rows(., nalx_prev_rx_opd_yr_tbl, ss_prev_rx_opd_yr_tbl,
            pg_prev_rx_opd_yr_tbl, allint_prev_rx_opd_yr_tbl) %>% 
    bind_rows(., prop_opioids_rx_target_tbl %>% 
              select(year, value) %>% 
              mutate(intervention = "Target"))

prev_rx_opd_mon_tbl$year_mon <- factor(prev_rx_opd_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

prop_opioids_rx_target_tbl$year_mon <- factor(prop_opioids_rx_target_tbl$year_mon,
                                              levels = paste(rep(v_yrs_prev,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"),
                                              labels = paste(rep(v_yrs_prev,
                                                                 each = 12),
                                                             month.abb[1:12],
                                                             sep = "_"))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% c("BPO_ACUTE", "BPO_CHRONIC", "BPO_CANCER",
                               "BPO_OTHER", "BPO_PALLIATIVE")) %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_rx_opd_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of prescription opioid use") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.5 Prevalence of illicit use

Code
v_states_illicit <- c("BI_ILLICIT", "BR_ILLICIT")
noint_prev_illicituse_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_illicit)$prop_uncalib_tbl
noint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
nalx_prev_illicituse_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = v_states_illicit)$prop_uncalib_tbl
nalx_prev_illicituse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
ss_prev_illicituse_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = v_states_illicit)$prop_uncalib_tbl
ss_prev_illicituse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
pg_prev_illicituse_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = v_states_illicit)$prop_uncalib_tbl
pg_prev_illicituse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean
allint_prev_illicituse_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = v_states_illicit)$prop_uncalib_tbl
allint_prev_illicituse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_illicit)$prop_uncalib_wei_mean


prev_illicituse_mon_tbl <- noint_prev_illicituse_mon_tbl %>% 
  bind_rows(., nalx_prev_illicituse_mon_tbl, ss_prev_illicituse_mon_tbl,
            pg_prev_illicituse_mon_tbl, allint_prev_illicituse_mon_tbl)

prev_illicituse_yr_tbl <- noint_prev_illicituse_yr_tbl %>% 
  bind_rows(., nalx_prev_illicituse_yr_tbl, ss_prev_illicituse_yr_tbl,
            pg_prev_illicituse_yr_tbl, allint_prev_illicituse_yr_tbl) 

prev_illicituse_mon_tbl$year_mon <- factor(prev_illicituse_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% v_states_illicit) %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_illicituse_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid use") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.6 Prevalence of misuse

Code
noint_prev_rxmisuse_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BPO_MISUSE")$prop_uncalib_tbl
noint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
nalx_prev_rxmisuse_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = "BPO_MISUSE")$prop_uncalib_tbl
nalx_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
ss_prev_rxmisuse_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = "BPO_MISUSE")$prop_uncalib_tbl
ss_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
pg_prev_rxmisuse_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = "BPO_MISUSE")$prop_uncalib_tbl
pg_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean
allint_prev_rxmisuse_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = "BPO_MISUSE")$prop_uncalib_tbl
allint_prev_rxmisuse_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = "BPO_MISUSE")$prop_uncalib_wei_mean


prev_rxmisuse_mon_tbl <- noint_prev_rxmisuse_mon_tbl %>% 
  bind_rows(., nalx_prev_rxmisuse_mon_tbl, ss_prev_rxmisuse_mon_tbl,
            pg_prev_rxmisuse_mon_tbl, allint_prev_rxmisuse_mon_tbl)

prev_rxmisuse_yr_tbl <- noint_prev_rxmisuse_yr_tbl %>% 
  bind_rows(., nalx_prev_rxmisuse_yr_tbl, ss_prev_rxmisuse_yr_tbl,
            pg_prev_rxmisuse_yr_tbl, allint_prev_rxmisuse_yr_tbl) 

prev_rxmisuse_mon_tbl$year_mon <- factor(prev_rxmisuse_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% "BPO_MISUSE") %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_rxmisuse_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid misuse") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.7 Prevalence of overdose

Code
v_states_od <- c("BO_OD_RX", "BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_od_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_od)$prop_uncalib_tbl
noint_prev_od_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
nalx_prev_od_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = v_states_od)$prop_uncalib_tbl
nalx_prev_od_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
ss_prev_od_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = v_states_od)$prop_uncalib_tbl
ss_prev_od_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
pg_prev_od_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = v_states_od)$prop_uncalib_tbl
pg_prev_od_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_od)$prop_uncalib_wei_mean
allint_prev_od_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = v_states_od)$prop_uncalib_tbl
allint_prev_od_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_od)$prop_uncalib_wei_mean


prev_od_mon_tbl <- noint_prev_od_mon_tbl %>% 
  bind_rows(., nalx_prev_od_mon_tbl, ss_prev_od_mon_tbl,
            pg_prev_od_mon_tbl, allint_prev_od_mon_tbl)

prev_od_yr_tbl <- noint_prev_od_yr_tbl %>% 
  bind_rows(., nalx_prev_od_yr_tbl, ss_prev_od_yr_tbl,
            pg_prev_od_yr_tbl, allint_prev_od_yr_tbl) 

prev_od_mon_tbl$year_mon <- factor(prev_od_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% v_states_od) %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_od_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of opioid-related overdose") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.7.1 prevalence of illicit overdose

Code
v_states_illicitod <- c("BO_OD_ILLICIT", "BR_OD_ILLICIT")
noint_prev_illicitod_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_illicitod)$prop_uncalib_tbl
noint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = v_states_illicitod)$prop_uncalib_wei_mean
nalx_prev_illicitod_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = v_states_illicitod)$prop_uncalib_tbl
nalx_prev_illicitod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = v_states_illicitod)$prop_uncalib_wei_mean
ss_prev_illicitod_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = v_states_illicitod)$prop_uncalib_tbl
ss_prev_illicitod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = v_states_illicitod)$prop_uncalib_wei_mean
pg_prev_illicitod_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = v_states_illicitod)$prop_uncalib_tbl
pg_prev_illicitod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = v_states_illicitod)$prop_uncalib_wei_mean
allint_prev_illicitod_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = v_states_illicitod)$prop_uncalib_tbl
allint_prev_illicitod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = v_states_illicitod)$prop_uncalib_wei_mean


prev_illicitod_mon_tbl <- noint_prev_illicitod_mon_tbl %>% 
  bind_rows(., nalx_prev_illicitod_mon_tbl, ss_prev_illicitod_mon_tbl,
            pg_prev_illicitod_mon_tbl, allint_prev_illicitod_mon_tbl)

prev_illicitod_yr_tbl <- noint_prev_illicitod_yr_tbl %>% 
  bind_rows(., nalx_prev_illicitod_yr_tbl, ss_prev_illicitod_yr_tbl,
            pg_prev_illicitod_yr_tbl, allint_prev_illicitod_yr_tbl) 

prev_illicitod_mon_tbl$year_mon <- factor(prev_illicitod_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% v_states_illicitod) %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_illicitod_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of illicit opioid-related overdose") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

10.6.7.2 prevalence of rx overdose

Code
noint_prev_rxod_mon_tbl <-  prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_OD_RX")$prop_uncalib_tbl
noint_prev_rxod_yr_tbl <- prev_fun(mod = mod_basecase, i = 1,
                                       v_states = "BO_OD_RX")$prop_uncalib_wei_mean
nalx_prev_rxod_mon_tbl <-  prev_fun(mod = mod_nalx_1, i = 2,
                                       v_states = "BO_OD_RX")$prop_uncalib_tbl
nalx_prev_rxod_yr_tbl <- prev_fun(mod = mod_nalx_1, i =2,
                                       v_states = "BO_OD_RX")$prop_uncalib_wei_mean
ss_prev_rxod_mon_tbl <-  prev_fun(mod = mod_ss_1, i = 3,
                                       v_states = "BO_OD_RX")$prop_uncalib_tbl
ss_prev_rxod_yr_tbl <- prev_fun(mod = mod_ss_1, i =3,
                                       v_states = "BO_OD_RX")$prop_uncalib_wei_mean
pg_prev_rxod_mon_tbl <-  prev_fun(mod = mod_pres_guid_1, i = 4,
                                       v_states = "BO_OD_RX")$prop_uncalib_tbl
pg_prev_rxod_yr_tbl <- prev_fun(mod = mod_pres_guid_1, i =4,
                                       v_states = "BO_OD_RX")$prop_uncalib_wei_mean
allint_prev_rxod_mon_tbl <-  prev_fun(mod = mod_all_int_1, i = 5,
                                       v_states = "BO_OD_RX")$prop_uncalib_tbl
allint_prev_rxod_yr_tbl <- prev_fun(mod = mod_all_int_1, i =5,
                                       v_states = "BO_OD_RX")$prop_uncalib_wei_mean


prev_rxod_mon_tbl <- noint_prev_rxod_mon_tbl %>% 
  bind_rows(., nalx_prev_rxod_mon_tbl, ss_prev_rxod_mon_tbl,
            pg_prev_rxod_mon_tbl, allint_prev_rxod_mon_tbl)

prev_rxod_yr_tbl <- noint_prev_rxod_yr_tbl %>% 
  bind_rows(., nalx_prev_rxod_yr_tbl, ss_prev_rxod_yr_tbl,
            pg_prev_rxod_yr_tbl, allint_prev_rxod_yr_tbl) 

prev_rxod_mon_tbl$year_mon <- factor(prev_rxod_mon_tbl$year_mon,
                                              levels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"
                                                  ), 5),
                                              labels = rep(
                                                paste(
                                                  rep(v_yrs_prev,
                                                      each = 12),
                                                  month.abb[1:12],
                                                  sep = "_"), 5))

ggplotly(ggplot(trace_tbl %>% 
         filter(state %in% "BO_OD_RX") %>%
           filter(scenario == "Scenario 1") %>% 
           group_by(cycle_num, mod) %>% 
           summarize(count = sum(count)) %>% 
           ungroup() %>% 
           left_join(., year_mon_cycle_tbl, by = c("cycle_num" = "cycle")),
         aes(x = cycle_num, y = count, colour = mod, key = year)) +
  geom_line() +
  ylab("Count") + xlab("Cycle Number") +
  scale_color_manual(values = colors_pal, name = "Intervention"))
Code
ggplotly(ggplot() +
  geom_point(data = prev_rxod_yr_tbl,
             aes(x = factor(year), y = value,
                 color = intervention, fill = intervention)) +
  xlab("Year") + ylab("Weighted Yearly Prevalence of rx opioid-related overdose") +
  scale_fill_discrete(name = "Intervention")+
  scale_color_discrete(name = "Intervention")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)))

11 Appendix 6: PSA results

Code
library(tidyverse)
library(plotly)
library(flextable)
theme_set(theme_minimal())
library(here)

m_outcomes_psa <- readRDS(here("01_data/m_outcomes_psa.RDS"))
point_est <- readRDS(here("01_data/point_estiamte.RDS"))
total_death_oddeaths <- readRDS(here("01_data/total_deaths_oddeaths.RDS"))


mod_names <- c("No Interventions", "Naloxone", "Safer Supply", "Prescription Guidelines", "All Interventions")
mod_names <- factor(mod_names, levels = mod_names, labels = mod_names)

11.1 Distribution of change in primary outcomes compared to No intervention

11.1.2 Overall deaths

Code
ggplotly(ggplot(data = diff_psa_dt %>% 
         filter(grp == "death"), aes(x = diff, fill = Intervention, color = Intervention)) + 
  geom_density(alpha = 0.9)+
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    xlab("Change in Deaths Compared to No Intervention") + ylab("") +
  geom_vline(data = point_est_diff %>%
                 filter(grp == "death"), aes(xintercept = diff, linetype = type))+
    labs(linetype = "") +
  facet_wrap(~Intervention, scales = "free_y") +
    theme_minimal())

11.1.3 Costs

Code
ggplotly(ggplot(data = diff_psa_dt %>% 
         filter(grp == "cost"), aes(x = diff, fill = Intervention, color = Intervention)) + 
  geom_density(alpha = 0.9)+
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    xlab("Change in Costs Compared to No Intervention (in millions)") + ylab("") +
  geom_vline(data = point_est_diff %>%
                 filter(grp == "cost"), aes(xintercept = diff/1000000, linetype = type), linewidth = 0.5)+
    labs(linetype = "") +
  facet_wrap(~Intervention, scales = "free_y") +
    theme_minimal())

11.2 Distribution of percent change in primary outcomes compared to No intervention

Code
diff_per_psa_dt <- outcomes_psa_dt %>% 
  mutate(cost_diff_per_nalox = cost_diff_per_nalox,
         cost_diff_per_ss = cost_diff_per_ss,
         cost_diff_per_pg = cost_diff_per_pg,
         cost_diff_per_all = cost_diff_per_all) %>% 
  select(cost_diff_per_nalox, cost_diff_per_ss, cost_diff_per_pg, cost_diff_per_all,
         death_diff_per_nalox, death_diff_per_ss, death_diff_per_pg, death_diff_per_all,
         oddeath_diff_per_nalox, oddeath_diff_per_ss, oddeath_diff_per_pg, oddeath_diff_per_all) %>% 
  pivot_longer(1:12, names_to = "group", values_to = "diff_per") %>% 
  separate(col = group, into = c("grp", "type1", "type2", "Intervention"), sep = "_") %>% 
  select(-type1, -type2)  %>% 
  mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
                               ifelse(Intervention == "ss", "Safer Supply",
                                      ifelse(Intervention == "pg", "Prescription Guidelines", 
                                             ifelse(Intervention == "all", "All Interventions", Intervention)))))


point_est_diff_per <- point_est %>% 
  select(cost_diff_per, deaths_diff_per, od_deaths_diff_per) %>% 
  cbind.data.frame(Intervention = mod_names[-1]) %>% 
  pivot_longer(1:3,names_to = "grp", values_to = "diff") %>% 
  separate(grp, into = c("grp", "type"), sep = "_") %>% select(-type) %>% 
  mutate(grp = ifelse(grp == "od", "oddeath", 
                      ifelse(grp == "deaths", "death", grp)),
         type = "Point Estimate")

11.2.2 Overall deaths

Code
ggplotly(ggplot(data = diff_per_psa_dt %>% 
         filter(grp == "death"), aes(x = diff_per, fill = Intervention, color = Intervention)) + 
  geom_density(alpha = 0.9)+
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    xlab("% Change in Deaths Compared to No Intervention") + ylab("") +
     geom_vline(data = point_est_diff_per %>%
                 filter(grp == "death"), aes(xintercept = diff, linetype = type))+
    labs(linetype = "") +
  facet_wrap(~Intervention, scales = "free_y") +
    theme_minimal())

11.2.3 Costs

Code
ggplotly(ggplot(data = diff_per_psa_dt %>% 
         filter(grp == "cost"), aes(x = diff_per, fill = Intervention, color = Intervention)) + 
  geom_density(alpha = 0.9)+
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    xlab("% Change in Costs Compared to No Intervention") + ylab("") +
   geom_vline(data = point_est_diff_per %>%
                 filter(grp == "cost"), aes(xintercept = diff, linetype = type))+
    labs(linetype = "") +
  facet_wrap(~Intervention, scales = "free_y") +
    theme_minimal())

11.3 Table

Code
mean_cost_diff_nalox <- mean(m_outcomes_psa[, "cost_diff_nalox"])
ci_cost_diff_nalox <- quantile(m_outcomes_psa[, "cost_diff_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_ss <- mean(m_outcomes_psa[, "cost_diff_ss"])
ci_cost_diff_ss <- quantile(m_outcomes_psa[, "cost_diff_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_pg <- mean(m_outcomes_psa[, "cost_diff_pg"])
ci_cost_diff_pg <- quantile(m_outcomes_psa[, "cost_diff_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_all <- mean(m_outcomes_psa[, "cost_diff_all"])
ci_cost_diff_all <- quantile(m_outcomes_psa[, "cost_diff_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_cost_diff <- c(paste0(round(ci_cost_diff_nalox[3]/1000000, 0), " [",
                           round(ci_cost_diff_nalox[1]/1000000, 0),", ",
                           round(ci_cost_diff_nalox[2]/1000000, 0), "]"),
                    paste0(round(ci_cost_diff_ss[3]/1000000, 0), " [",
                           round(ci_cost_diff_ss[1]/1000000, 0),", ",
                           round(ci_cost_diff_ss[2]/1000000, 0), "]"),
                    paste0(round(ci_cost_diff_pg[3]/1000000, 0), " [",
                           round(ci_cost_diff_pg[1]/1000000, 0),", ",
                           round(ci_cost_diff_pg[2]/1000000, 0), "]"),
                    paste0(round(ci_cost_diff_all[3]/1000000, 0), " [",
                           round(ci_cost_diff_all[1]/1000000, 0),", ",
                           round(ci_cost_diff_all[2]/1000000, 0), "]"))


mean_death_diff_nalox <- mean(m_outcomes_psa[, "death_diff_nalox"])
ci_death_diff_nalox <- quantile(m_outcomes_psa[, "death_diff_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_ss <- mean(m_outcomes_psa[, "death_diff_ss"])
ci_death_diff_ss <- quantile(m_outcomes_psa[, "death_diff_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_pg <- mean(m_outcomes_psa[, "death_diff_pg"])
ci_death_diff_pg <- quantile(m_outcomes_psa[, "death_diff_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_all <- mean(m_outcomes_psa[, "death_diff_all"])
ci_death_diff_all <- quantile(m_outcomes_psa[, "death_diff_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_death_diff <- c(paste0(round(ci_death_diff_nalox[3], 0), " [",
                           round(ci_death_diff_nalox[1], 0),", ",
                           round(ci_death_diff_nalox[2], 0), "]"),
                    paste0(round(ci_death_diff_ss[3], 0), " [",
                           round(ci_death_diff_ss[1], 0),", ",
                           round(ci_death_diff_ss[2], 0), "]"),
                    paste0(round(ci_death_diff_pg[3], 0), " [",
                           round(ci_death_diff_pg[1], 0),", ",
                           round(ci_death_diff_pg[2], 0), "]"),
                    paste0(round(ci_death_diff_all[3], 0), " [",
                           round(ci_death_diff_all[1], 0),", ",
                           round(ci_death_diff_all[2], 0), "]"))

mean_oddeath_diff_nalox <- mean(m_outcomes_psa[, "oddeath_diff_nalox"])
ci_oddeath_diff_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_ss <- mean(m_outcomes_psa[, "oddeath_diff_ss"])
ci_oddeath_diff_ss <- quantile(m_outcomes_psa[, "oddeath_diff_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_pg <- mean(m_outcomes_psa[, "oddeath_diff_pg"])
ci_oddeath_diff_pg <- quantile(m_outcomes_psa[, "oddeath_diff_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_all <- mean(m_outcomes_psa[, "oddeath_diff_all"])
ci_oddeath_diff_all <- quantile(m_outcomes_psa[, "oddeath_diff_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_oddeath_diff <- c(paste0(round(ci_oddeath_diff_nalox[3], 0), " [",
                           round(ci_oddeath_diff_nalox[1], 0),", ",
                           round(ci_oddeath_diff_nalox[2], 0), "]"),
                    paste0(round(ci_oddeath_diff_ss[3], 0), " [",
                           round(ci_oddeath_diff_ss[1], 0),", ",
                           round(ci_oddeath_diff_ss[2], 0), "]"),
                    paste0(round(ci_oddeath_diff_pg[3], 0), " [",
                           round(ci_oddeath_diff_pg[1], 0),", ",
                           round(ci_oddeath_diff_pg[2], 0), "]"),
                    paste0(round(ci_oddeath_diff_all[3], 0), " [",
                           round(ci_oddeath_diff_all[1], 0),", ",
                           round(ci_oddeath_diff_all[2], 0), "]"))


effects_tbl <- cbind.data.frame(mod_names[-1], mean_cost_diff,
                                mean_death_diff, mean_oddeath_diff)

effects_tbl |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         mean_cost_diff = "Discounted Net Present \n Costs in Millions ",
                         mean_death_diff = "Deaths ",
                         mean_oddeath_diff = "Opioid-related Overdose Deaths")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Change Compared to No Intervention
Median [95% Credible Intervals]

Interventions

Discounted Net Present
Costs in Millions

Deaths

Opioid-related Overdose Deaths

Naloxone

209 [107, 362]

625 [414, 926]

-10089 [-15009, -6697]

Safer Supply

4216 [3651, 4821]

-14180 [-17532, -11074]

-3221 [-4105, -2529]

Prescription Guidelines

-25455 [-31480, -19866]

14417 [8610, 20565]

-5458 [-11045, -1624]

All Interventions

-21058 [-27059, -15403]

960 [-5687, 7753]

-18544 [-28496, -11309]

Code
mean_cost_diff_per_nalox <- mean(m_outcomes_psa[, "cost_diff_per_nalox"])
ci_cost_diff_per_nalox <- quantile(m_outcomes_psa[, "cost_diff_per_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_per_ss <- mean(m_outcomes_psa[, "cost_diff_per_ss"])
ci_cost_diff_per_ss <- quantile(m_outcomes_psa[, "cost_diff_per_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_per_pg <- mean(m_outcomes_psa[, "cost_diff_per_pg"])
ci_cost_diff_per_pg <- quantile(m_outcomes_psa[, "cost_diff_per_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_cost_diff_per_all <- mean(m_outcomes_psa[, "cost_diff_per_all"])
ci_cost_diff_per_all <- quantile(m_outcomes_psa[, "cost_diff_per_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_cost_diff_per <- c(paste0(round(ci_cost_diff_per_nalox[3], 2), " [",
                           round(ci_cost_diff_per_nalox[1], 2),", ",
                           round(ci_cost_diff_per_nalox[2], 2), "]%"),
                    paste0(round(ci_cost_diff_per_ss[3], 2), " [",
                           round(ci_cost_diff_per_ss[1], 2),", ",
                           round(ci_cost_diff_per_ss[2], 2), "]%"),
                    paste0(round(ci_cost_diff_per_pg[3], 2), " [",
                           round(ci_cost_diff_per_pg[1], 2),", ",
                           round(ci_cost_diff_per_pg[2], 2), "]%"),
                    paste0(round(ci_cost_diff_per_all[3], 2), " [",
                           round(ci_cost_diff_per_all[1], 2),", ",
                           round(ci_cost_diff_per_all[2], 2), "]%"))


mean_death_diff_per_nalox <- mean(m_outcomes_psa[, "death_diff_per_nalox"])
ci_death_diff_per_nalox <- quantile(m_outcomes_psa[, "death_diff_per_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_per_ss <- mean(m_outcomes_psa[, "death_diff_per_ss"])
ci_death_diff_per_ss <- quantile(m_outcomes_psa[, "death_diff_per_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_per_pg <- mean(m_outcomes_psa[, "death_diff_per_pg"])
ci_death_diff_per_pg <- quantile(m_outcomes_psa[, "death_diff_per_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_death_diff_per_all <- mean(m_outcomes_psa[, "death_diff_per_all"])
ci_death_diff_per_all <- quantile(m_outcomes_psa[, "death_diff_per_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_death_diff_per <- c(paste0(round(ci_death_diff_per_nalox[3], 2), " [",
                           round(ci_death_diff_per_nalox[1], 2),", ",
                           round(ci_death_diff_per_nalox[2], 2), "]%"),
                    paste0(round(ci_death_diff_per_ss[3], 2), " [",
                           round(ci_death_diff_per_ss[1], 2),", ",
                           round(ci_death_diff_per_ss[2], 2), "]%"),
                    paste0(round(ci_death_diff_per_pg[3], 2), " [",
                           round(ci_death_diff_per_pg[1], 2),", ",
                           round(ci_death_diff_per_pg[2], 2), "]%"),
                    paste0(round(ci_death_diff_per_all[3], 2), " [",
                           round(ci_death_diff_per_all[1], 2),", ",
                           round(ci_death_diff_per_all[2], 2), "]%"))

mean_oddeath_diff_per_nalox <- mean(m_outcomes_psa[, "oddeath_diff_per_nalox"])
ci_oddeath_diff_per_nalox <- quantile(m_outcomes_psa[, "oddeath_diff_per_nalox"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_per_ss <- mean(m_outcomes_psa[, "oddeath_diff_per_ss"])
ci_oddeath_diff_per_ss <- quantile(m_outcomes_psa[, "oddeath_diff_per_ss"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_per_pg <- mean(m_outcomes_psa[, "oddeath_diff_per_pg"])
ci_oddeath_diff_per_pg <- quantile(m_outcomes_psa[, "oddeath_diff_per_pg"],
                                probs = c(0.025, 0.975, 0.5))

mean_oddeath_diff_per_all <- mean(m_outcomes_psa[, "oddeath_diff_per_all"])
ci_oddeath_diff_per_all <- quantile(m_outcomes_psa[, "oddeath_diff_per_all"],
                                probs = c(0.025, 0.975, 0.5))


mean_oddeath_diff_per <- c(paste0(round(ci_oddeath_diff_per_nalox[3], 2), " [",
                           round(ci_oddeath_diff_per_nalox[1], 2),", ",
                           round(ci_oddeath_diff_per_nalox[2], 2), "]%"),
                    paste0(round(ci_oddeath_diff_per_ss[3], 2), " [",
                           round(ci_oddeath_diff_per_ss[1], 2),", ",
                           round(ci_oddeath_diff_per_ss[2], 2), "]%"),
                    paste0(round(ci_oddeath_diff_per_pg[3], 2), " [",
                           round(ci_oddeath_diff_per_pg[1], 2),", ",
                           round(ci_oddeath_diff_per_pg[2], 2), "]%"),
                    paste0(round(ci_oddeath_diff_per_all[3], 2), " [",
                           round(ci_oddeath_diff_per_all[1], 2),", ",
                           round(ci_oddeath_diff_per_all[2], 2), "]%"))


effects_tbl_per <- cbind.data.frame(mod_names[-1], mean_cost_diff_per,
                                mean_death_diff_per, mean_oddeath_diff_per)

effects_tbl_per |> 
  flextable() |> 
  add_header_row(colwidths = c(1, 3),
                      values = c("", "Percent Change Compared to No Intervention \n Median [95% Credible Intervals]")) |>
  set_header_labels(values = list(
                         `mod_names[-1]` = "Interventions",
                         mean_cost_diff_per = "Discounted Net Present Costs",
                         mean_death_diff_per = "Deaths ",
                         mean_oddeath_diff_per = "Opioid-related Overdose Deaths")) |>
  theme_vanilla() |>
  set_table_properties(layout = "autofit") |>
  align_text_col(align = "center", header = TRUE, footer = TRUE)

Percent Change Compared to No Intervention
Median [95% Credible Intervals]

Interventions

Discounted Net Present Costs

Deaths

Opioid-related Overdose Deaths

Naloxone

0.03 [0.01, 0.05]%

0.01 [0.01, 0.02]%

-3.46 [-4.13, -3.13]%

Safer Supply

0.54 [0.46, 0.61]%

-0.31 [-0.38, -0.24]%

-1.11 [-1.97, -0.66]%

Prescription Guidelines

-3.23 [-3.94, -2.55]%

0.32 [0.19, 0.45]%

-1.81 [-2.76, -0.93]%

All Interventions

-2.67 [-3.39, -1.98]%

0.02 [-0.13, 0.17]%

-6.43 [-7.16, -5.58]%

Code
tbl_effect_plot1 <- cbind.data.frame(value = c("lb", "ub", "median"),
                               ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
                               ci_cost_diff_per_pg, ci_cost_diff_per_all,
                               ci_death_diff_per_nalox, ci_death_diff_per_ss,
                               ci_death_diff_per_pg, ci_death_diff_per_all,
                               ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
                               ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>% 
  pivot_longer(2:13, names_to = "grp", values_to = "per_diff") %>% 
  pivot_wider(names_from = "value", values_from = "per_diff") %>% 
  separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>% 
  select(-c(type1, type2, type3)) %>% 
  mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
                               ifelse(Intervention == "ss", "Safer Supply",
                                      ifelse(Intervention == "pg", "Prescription Guidelines", 
                                             ifelse(Intervention == "all", "All Interventions", Intervention)))))
  
tbl_effect_plot1$Intervention <- factor(tbl_effect_plot1$Intervention, levels = mod_names, labels = mod_names)

ggplotly(ggplot(data = tbl_effect_plot1, aes(x = group, y = median))+
  geom_jitter(aes(color = Intervention, shape = Intervention), size = 2, position = position_dodge(0.25)) +
    # geom_errorbar(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), width = 0.3, position = position_dodge(0.1)) +
    geom_pointrange(aes(ymin = lb, ymax = ub, color = Intervention, shape = Intervention), position = position_dodge(0.25)) + 
  scale_color_brewer(palette = "Dark2") +
  scale_shape_manual(values = c(15, 16, 17, 18)) +
  # scale_shape_manual(values = c(15, 16)) +
    geom_hline(yintercept = 0) +
  ylab("Percent Change Compared to No Intervention (%)") +
  xlab("") +
  scale_x_discrete(name ="", 
                    labels=c("Costs","Deaths","Opioid-related Overdose Deaths"))+
    theme_minimal())
Code
tbl_effect_plot2 <- cbind.data.frame(value = c("cost_lb", "cost_ub", "cost_median"),
                                     ci_cost_diff_per_nalox, ci_cost_diff_per_ss,
                                     ci_cost_diff_per_pg, ci_cost_diff_per_all) %>% 
  pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
   pivot_wider(names_from = "value", values_from = "per_diff") %>% 
  separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>% 
  select(-c(type1, type2, type3, group)) %>% 
  full_join(., cbind.data.frame(value = c("oddeath_lb", "oddeath_ub", "oddeath_median"),
                                     ci_oddeath_diff_per_nalox, ci_oddeath_diff_per_ss,
                                     ci_oddeath_diff_per_pg, ci_oddeath_diff_per_all) %>% 
              pivot_longer(2:5, names_to = "grp", values_to = "per_diff") %>%
              pivot_wider(names_from = "value", values_from = "per_diff") %>% 
              separate(grp, into = c("type1", "group", "type2", "type3", "Intervention")) %>% 
              select(-c(type1, type2, type3, group)), by = "Intervention") %>% 
  mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
                               ifelse(Intervention == "ss", "Safer Supply",
                                      ifelse(Intervention == "pg", "Prescription Guidelines", 
                                             ifelse(Intervention == "all", "All Interventions", Intervention)))))
  
  tbl_effect_plot2$Intervention <- factor(tbl_effect_plot2$Intervention, levels = mod_names, labels = mod_names)


ggplotly(ggplot(data = tbl_effect_plot2, aes(x = cost_median, y = oddeath_median)) +
  geom_point(aes(color = Intervention), size = 1) +
    geom_pointrange(aes(ymin = oddeath_lb, ymax = oddeath_ub, color = Intervention)) +
  geom_errorbarh(aes(xmax = cost_lb, xmin = cost_ub,color = Intervention),  height = 0) +
  scale_color_brewer(palette = "Dark2") +
  # scale_shape_manual(values = c(15, 16, 17, 18)) +
    geom_hline(yintercept = 0, linewidth = 0.25) +
  geom_vline(xintercept = 0, linewidth = 0.25) +
  xlab("Change in Costs Compared to No Intervention (%)") +
  ylab("Change in Opioid-related Overdose Deaths \n Compared to No Intervention (%)") +
  theme_minimal())

11.4 Posterior distributions for total Deaths and total OD-deaths

Code
death_psa_dt <- outcomes_psa_dt %>% 
  select(tot_death_no, tot_death_nalox, tot_death_ss, tot_death_pg, tot_death_all,
         tot_oddeath_no, tot_oddeath_nalox, tot_oddeath_ss, tot_oddeath_pg, tot_oddeath_all) %>% 
  pivot_longer(1:10, names_to = "group", values_to = "total") %>% 
  separate(col = group, into = c("type", "grp", "Intervention"), sep = "_") %>% 
  select(-type) %>% 
  mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
                               ifelse(Intervention == "ss", "Safer Supply",
                                      ifelse(Intervention == "pg", "Prescription Guidelines", 
                                             ifelse(Intervention == "all", "All Interventions",
                                                    ifelse(Intervention == "no", "No Interventions", Intervention))))))
death_psa_dt$Intervention <- factor(death_psa_dt$Intervention,
                                    levels = mod_names, labels = mod_names)

total_death_oddeaths <- total_death_oddeaths %>% 
  pivot_longer(3:7, names_to = "Intervention", values_to = "total") %>% 
  mutate(grp = ifelse(state == "BO_OD_DEATH", "oddeath", "death"),
         type = "Point Estimate") %>% 
  select(-state)

ggplot(data = death_psa_dt %>% 
         filter(grp == "oddeath"), aes(x = total, fill = Intervention, color = Intervention)) + 
  geom_density(alpha = 0.9)+
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    xlab("Total Opioid-related Overdose Deaths after 15 years") + ylab("") +
    geom_vline(data = total_death_oddeaths %>%
                 filter(grp == "oddeath"), aes(xintercept = total, linetype = type), linewidth = 0.5)+
    labs(linetype = "") +
  facet_wrap(~Intervention, scales = "free_y") +
    theme_minimal()

Code
ggplot(data = death_psa_dt %>% 
         filter(grp == "death"), aes(x = total/1000, fill = Intervention, color = Intervention)) + 
  geom_density(alpha = 0.9)+
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2") +
    xlab("Total Deaths after 15 years (in thousands)") + ylab("") +
  geom_vline(data = total_death_oddeaths %>%
                 filter(grp == "death"), aes(xintercept = total/1000, linetype = type), linewidth = 0.5)+
    labs(linetype = "") +
  facet_wrap(~Intervention, scales = "free_y") +
    theme_minimal()

11.5 New OD Deaths per year

Code
outcomes_psa_dt_inci <- outcomes_psa_dt %>% 
  select(35:104) 



inci_od_dt <- data.frame(nm = rep(NA, length(outcomes_psa_dt_inci)),
           med = rep(NA, length(outcomes_psa_dt_inci)),
           lb = rep(NA, length(outcomes_psa_dt_inci)),
           ub = rep(NA, length(outcomes_psa_dt_inci)))


for (i in 1:length(outcomes_psa_dt_inci)){
  inci_od_dt$nm[i] <- names(outcomes_psa_dt_inci)[i]
  inci_od_dt[i, 2:4] <- as.numeric(quantile(outcomes_psa_dt_inci[i, ], c(0.5, 0.025, 0.975)))
  
}


inc_od_death_tbl_mod <- readRDS(here("01_data/inc_od_death_tbl_mod.RDS")) %>% 
  filter(scenario == "Scenario 1") %>% 
  select(-scenario) %>% 
  rename(Intervention = intervention,
         Year = year, 
         value = inci_od_deaths) %>% 
  mutate(grp = "Point Estimate")

inci_od_dt_plot <- inci_od_dt %>% 
  separate(nm, into = c("type", "type1", "Intervention", "Year"), sep = "_") %>% 
  select(-c(type, type1)) %>%
  rename(value = med) %>% 
  mutate(Intervention = ifelse(Intervention == "nalox", "Naloxone",
                               ifelse(Intervention == "ss", "Safer Supply",
                                      ifelse(Intervention == "pg", "Prescription Guidelines", 
                                             ifelse(Intervention == "all", "All Interventions",
                                                    ifelse(Intervention == "no", "No Interventions", Intervention))))),
         Year = c(rep(2016:2029, 5)),
         grp = "PSA results") %>% 
  bind_rows(., inc_od_death_tbl_mod) %>% 
  mutate(newgrp = paste0(Intervention, "_", grp))

inci_od_dt_plot$Intervention <- factor(inci_od_dt_plot$Intervention, levels = mod_names, labels = mod_names)


inc_od_death_tbl_target <- readRDS(here("01_data/inc_od_death_tbl_target.RDS")) %>% 
  mutate(grp = "Target")




ggplotly(ggplot(data = inci_od_dt_plot, aes(x = Year, y = value, color = grp, fill = grp)) +
  geom_line() +
  geom_point(data = inc_od_death_tbl_target, aes(x = year, y = inci_od_deaths, color = grp, fill = grp))+
  scale_color_manual(values = RColorBrewer::brewer.pal(3, "Dark2")[c(1,3,2)]) +
    scale_fill_manual(values = RColorBrewer::brewer.pal(3, "Dark2")[c(1,3, 2)]) +
    # scale_color_brewer(palette = "Dark2") +
  # scale_fill_brewer(palette = "Dark2") +
    geom_ribbon(aes(ymin = lb, ymax = ub), alpha = 0.5) +
  facet_wrap(~Intervention, scales = "free") +
    theme(legend.title = element_blank()))